In the name of Allah, 
Most Gracious, Most Merciful. 




Article Title: 

Miniatoric Infinity Plot And lt,s Applications 




Digest Description Of Miniatoric Infinity Plot Concept : 

• If either of the range end points of the horizontal 
range contains ±infinity, an infinity plot is generated. 

• An infinity plot is obtained by transforming -infinity .. 
-i-infinity to \ by a transformation that 
approximates arctan. This is a nice way of getting the 
entire picture of f(x) on the display. 

• Such a graph, although distorted near x = -infinity 
and + infinity, contains a lot of information about the 
features of f(x). 
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In the name of Allah, 



Most Gracious, Most Merciful. 






Blessed Martyrdom Of My Beautiful And Kind 
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In the name of Allah, 

Most Gracious, Most Merciful. 






The Sky Book ( Quran ) : Surat ( Al-Ahzab , or The Confederates ) : Verse (23) : 

^ ^ 

23 < |ii ^LjJo j 


Among the Believers are men who have been true to their covenant with 
Allah: of them some have completed their vow (to the extreme), and 
some (still) wait: but they have never changed (their determination) in 
the least: 




The Sky Book ( Quran ) : Surat ( Al-Baqarah ; or The Heifer ) : Verse (154) : 

4 154 < o^aJ) T ( jSsi j jJj dS4\ ^JwcLAj ^ J 24 V 


And say not of those who are slain in the way of Allah: "They are dead." 
Nay, they are living, though ye perceive (it) not. 


Introduction 


T ry finding the limit of the following function using Mathematica Software , which is 

known as one of the most powerful CAS "Computer Algebra Systems" when x — > +oo. 

[xl 

f(x) = — , [x] = intg(x) ( [x] = n <h> n<x<n + l,n6Z) 

X 

lntg=Function[ {x}, lf[ lntegerPart[x]>=0 , lntegerPart[x], lntegerPart[x]-l ]]; 

You will observe that Mathematica is unable to solve the limit . 

Attention : 

News from new versions of mathematica, year 2017: 

I sent my discovery to all CAS " Computer Algebra Systems " That I knew . 

Mathematica is now able to solve the limit after my discovery . 


In this article we will develop methods to find such limits , one method is graphical And the 
other method will be analytical . The graphical method is plotting the transformed function of 
F[x] in the domain — ^<x<+^as xis varying in the range -oo < x < +oo . We name such 

plots as Infinity Plots . Such Plots are introduced in Maple Software But there is no extended 
explanation how we can do this .( I searched the internet to find a solution and finally I found a 
method ) . 

Maple Software description about Infinity Plots : 


Calling Sequence 

plot(f, h, v, options) 

Parameters 


f(x) - 

acceptable function 

h 

horizontal range 

V 

vertical range (optional) 

Description 
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If either of the range end points of the horizontal range contains +-infinity, an 
infinity plot is generated. 

An infinity plot is obtained by transforming -infinity .. infinity to — ^ < x < + ^ 

by a transformation that approximates ArcTan. This is a nice way of getting the entire 
picture of f(x) on the display. Such a graph, although distorted near x = -infinity and 
infinity, contains a lot of information about the features of f(x). 

Because the view is already determined for infinity plots, the view option has no 

effect. 

Example : 


Mathematica Statements To Plot An Infinity Plot : 


F — Function 



G — Function 


{x}, ArcTan |V[Tan[x]]]J 


G[x] 

Plot[ArcTan[Cot[x]],{x,-(|),| }] 


The Result : 


Function 


w 4 


Fun cti on [{x}, ArcTan [F [Tan [x] ] ] ] 
ArcTan[Cot[x]] 


Plot[ArcTan[Cot[x]], {x, 


There are also other situations where these methods is helpful ,they are as follows : 

• Comparing the speed of different functions in growing to infinity . 

• Finding the finite number of roots of a function graphically as x is varying at the range 

—OO < X < +00 . 

• Plotting the asymptote of functions . 

• Developing some codes to find the finite number of roots of functions using the 
Kronecher - Picard Theory . 

• Developing a method to find the finite number of infinities of the function . 

• Developing a method to find the Absolute maximum and minimum of functions in the 
domain (—00, +00) . 

• 3 D Infinity plots 

• Developing some codes to find the solutions of system of nonlinear equations as the 
variables are changing in the range of -00 < x < +00 . 


Theory of infinity plots 



W hen plotting a function if either the range end points of horizontal range contains 
± infinity , an Infinity plot is generated. 

Suppose that we have the function f(x) and we want to infinity plot it . We can use 
the following transformation to get the entire picture of f(x) on the display . 

The transformation is : 


ArcTan (/('Tan(x))) 


D f = 




n n 
2 9 2 . 


If this transformation is plotted as x vary at the range — j < x < + 1 the f(x ) will be depicted 
on the display for x varying at the range — oo < x < +oo . 

Infinity plots are useful to inspect the behavior of one or more functions as — > ±oo . Scaling in 
both axes is not linear ; if the function to be plotted is oscillating for high |x|-values , the details 
are not shown correctly . In other words , Infinity plot is not suitable for plotting alternative 
functions. Note that it doesn,t make any sense to plot "Infiniy" plots together with normal plots 
in the same graph , since scaling is linear in normal plots , while this is not true in "Infinite" plots. 
Therefore we can not use "Infinity" plots together with the normal plots in the same graph. 

Here is some examples generated using Mathematica Software . To produce the plots two 
Mathematica statements are used , "Function" statement and "Plot" statement. 


Function Statement M : 

• Function [body] or body& is a pure function. The formal parameters are # (or #1), #2, etc. 

• Function [ x , body ] is a pure function with a single formal parameter x. 

• Function [ {xi, x 2 , ... }, body ] is a pure function with a list of formal parameters. 


Plot Statement M : 

• Plot [f, { x , x mm , x max } ] generates a plot of /as a function of x from x niin to x max . 

• Plot [ {/i, f 2 , ... } , {x, x mm , x max }] plots several functions /i. 
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Example 1 : 

fix) = x 2 

Mathematica Code : 


F= Function [{x},x A 2]; 

Plot[ArcTan[ F[Tan[x]]], {x, -n/2,n/2}]; 



Example 2 : 

f{x) = 5in(x) 

Mathematica Code : 


F=Function [{x},Sin[x ]]; 

Plot[ArcTan[ F[Tan[x]]], {x, -n/2,n/2}]; 
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Mathematica Code : 


F=Function [{x},5/(x A 2-4) ]; 
Plot[ArcTan[ F[Tan[x]]], {x, -n/2,n/2}]; 



Example 4 : 

f(x) = Ln(x) 
g(x) = e x 

Mathematica Code : 


F= Function [{x},Log[x]]; 

G=Function[ {x}, e x ] 

Plot[{ArcTan[F[Tan[x]]],ArcTan[G[Tan[x]]]},{x,-n:/2,n/2}]; 
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Infinity Plot Applications 


Limits Of Functions 


5 everal conclusions concerning the limits of the functions can be obtained using "Infinity" 
plots. Even there are cases where Mathematica which is known as the most powerful CAS 
"Computer Algebra System" , can not handle, but you can use this method to find the 
limits of functions graphically so easily . 

For Example if we try to find the limit of the following function as x -* +co using the 
Mathematica Software there will be no result . 

[x] 

f (x) = — 

x 

First we should define the function intg[x] for Mathematica : 


lntg=Function[ {x}, lf[ lntegerPart[x]>=0 , lntegerPart[x], lntegerPart[x]-l ]]; 


To define the function Intg we have used IntegerPart Statement which is a build in 
Mathematica function , which gives the Integer Part ofx, where x is the input variable. 

Now we try to find the limit of the defined function "Intg" divided byxasx—> +oo . 


Limit [ In tg[x]/x,x-> «> ] 


The output is : 


Lim[ 


If[IntegerPart[x] > 0 , IntegerPart[x], IntegerPart[x] — 1] 

x 


,x -* oo] 


As you can see there will be no result and Mathematica was unable to handle that limit. 
Now we are going to try the graphical method using the "Infinity" plots. 

Attention : 

News from new versions of mathematica, year 2017: 

I sent my discovery to all CAS " Computer Algebra Systems " That I knew . 

Mathematica is now able to solve the limit after my discovery . 



From the Mathematics courses we know that [x] as x -» +oo is equal to x .So we can guess that 
the limit will be simplified to : 



^ [x] 

To show that the above assumption is true we will plot f(x) = — and g(x) = 1 , on the same 
Infinity graph. 


Mathematica Code: 


lntg=Function[{x},lf[lntegerPart[x}>=0 ,lntegerPart[x],lntegerPart[x]-l]]; 
F=Function[ {xj, In tg[x ]/x]; 

Pl=Plot[{ArcTan[F[Tan[x]]]},{x,- tz /2, tt / 2},P!otRange->{- n /2, n/2}]; 
P2=Plot[{ArcTan[l]},{x,- n /2, n /2},P!otRange-*{- n /2, n/2}]; 
P3=Plot[{ArcTan[F[Tan[x]]],ArcTan[l]},{x,- n/2, n /2},PlotRange^>{- n/2, n/2} 
,PlotStyle^>{Thickness[.01 ],Dashing[{.01 }]}]; 



1.5 

1 


1.5 

1 


1.5 

1 


—\ 

0.5 


\ J S| ^ 0.5 


0.5 



- 1.5 -1 

- 0.5 

- 0.5 

-1 

- 1.5 

; o.5 

1 1.5 - 1.5 -1 - 0.5 

- 0.5 

-1 

- 1.5 

\ — 1 

LO 

\ — 1 

LO 

\ — 1 

LO 

o 



- 0.5 

- 0.5 

-1 

- 1.5 

; o.5 

1 1.5 


[ x ] 

From the graph you can see that — is getting near to 1 as x —> +oo . So the above Conclusion 
was true. 


Further Examples 
Example 1 : 

lim x ^ ±0O = ? , lim x ^_ 2 - 7Z ^ = ? , lim x ^_ 2+ 7Z ^ = ?, lim x ^ +2 - — =? . 


lim 


(x 2 - 4 ) 

X-» + 2+ f v 2_. 


(x 2 -4) 


(x 2 -4) 


(x 2 -4) 


= 7 


(x 2 — 4) 


Mathematica Code : 


F=Function }{x},5/(x A 2-4) ]; 

Plot [ {ArcTan[F[Tan[x}]]}, {x, -n/2, n/2}]; 




You con see that the results are : 
lim x 

lim 

lim 

lim 


5 

= 0 

lx ^±0° (X 2_ 4) 

5 

= +00 

lx ^- 2 " ( x 2 - 4 ) 

5 

= —00 

‘ x ^“ 2+ (x 2 — 4) 

5 

= —00 

l X^ + 2_ ( X 2_4) 

5 

= +00 

‘ x ^ + 2+ ( X 2_4) 


Example 2 : 

lim x ^ 0 + Ln(x) = ? 

Mathematica Code : 


F= Function [{x},Log[x]J; 
Plot[{ArcTan[F[Tan[x]]]},{x,-n/2,n/2}]; 
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You can see that the result is — oo . 


Example 3: 

l im x->+00 = 


? 


Mathematica Code : 


X 5 

F=Function[{x},— ]; 

Plot[ {ArcTan[F[Tan[x]]]}, {x, -n/2, n/2}] 



The result is 0. 

Checking the result using Mathematica : 


Limit[F[x],x-> Infinity] 


The result is 0 . 
Example 4: 

2 

lim x _, +00 (x + e x ) x =? 
Mathematica Code : 


F=Function[{x},(x + e x ) *]; 

Pl=Plot[ {ArcTan[ F[Tan[x]]]}, {x, -n/2, n /2], PlotRange^{1.4, n /2}]; 
P2=Plot[{ArcTan[e 2 ], {x,-n/2, n /2},P!otRange->{1.4, n/2}]; 


P3=Plot[{ArcTan[F[Tan[x]]],ArcTan[e 2 ], },{x,- n/2, n /2},PlotRange->{1.4, n/2}, 
PlotStyle^>{Thickness[.001 ],Dashing[{.01 }]}]; 

Show[GraphicsArray[{Pl,P2,P3}]]; 
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As you an see from the code and the graph the result is e 2 . 
Checking the result : 


Limit[ F[x],x->lnfinity] 


The result is e 2 . 

Example 5: 

lim x _, +00 (3x 2 + 4x)/(4x 2 + 1) =? 
Mathematica Code : 


F=Function[{x}, (3x 2 + 4x)/(4x 2 + 1)7; 

Pl=Plot[{ArcTan[F[Tan[x]]]},{x,- n /2, n /2},PlotRange^>{- n/2, n/2}]; 
P2=Plot[{ArcTan[3/4]},{x,- n /2, n /2},PiotRange->{- n /2, n /2},PlotStyle->Dashing[{.01}]]; 
P3=Plot[{ArcTan[F[Tan[x]]],ArcTan[3/4]},{x,- n/2, n/2},PlotRange^>{- n/2, n/2}, 
PlotStyle^>{Thickness[.001 ],Dashing[{.01 }]}]; 

Showf GraphicsArray[{Pl, P2, P3}]]; 

Limit [ F[x],x^>-lnfinity] 
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As you an see from the code and the graph the result is - . 

4 

Checking the result : 


Limit[F[x],x->lnfinity] 


3 

The result is 

4 


Example 6: 

x 


x->+oo Ln(x) x Ln(Ln(x)) 


Mathematica Code : 


Plot[{ArcTan[F[Tan[x]]]},{x,- n /2, n /2},P!otRange->{-n/2, n/2 }] 



The result is + oo . 


Example 7 : 


lim 

X— >+co 


Ln(x) 

yj—l + Lr^x)^/! + Ln(x) 


Mathematica Code : 


F=Function[{x},- 


Log(x) 


=]; 


’ N /-l+L°gfx) N /l+Log(x) 

Pl=Plot[{ArcTan[F[Tan[x]]]},{x,- n /2, n /2},PlotRange->{- n /2, n/2}]; 
P2=Plot[{ArcTan[l]},{x / - n /2, n /2},PlotRange^>{- n/2, n /2},PlotStyle^>Dashing[{.01}]]; 


P3=Plot[{ArcTan[F[Tan[x]]],ArcTan[l]},{x,- n/2, n /2},PlotRange^>{- n/2, n/2}, 
PlotStyle^>{Thickness[.001 ],Dashing[{.01 }]}]; 

Show[ GraphicsArray[{Pl, P2, P3}]]; 

Limit[ F[x],x^-lnfinity] 



As you an see from the code and the graph the result is 1 . 


Example 8 : 

3 X - 3 _x 

lim QX i Q— X 

x^+oo 3 X + 3 x 

Mathematica Code : 


F=Function[ M, 3X+3 _ X 7 ; 

Pl=Plot[{ArcTan[F[Tan[x]]]},{x,- n/2, n /2},PlotRange^>{- n/2, n/2}]; 
P2=Plot[{ArcTan[-l]},{x,- n/2, n /2},PlotRange^>{- n/2, n /2},PlotStyle^>Dashing[{.01}]]; 
P3=Plot[{ArcTan[F[Tan[x]]],ArcTan[-l]},{x,- n/2, n /2},P!otRange->{- n /2, n/2}, 
PlotStyle—>{Thickness[.001 ],Dashing[{.01 }]}]; 

Show [ GraphicsArray[{Pl, P2, P3}]]; 

Limit [ F[x],x^>-lnfinity] 



As you an see from the code and the graph the result is — 1 . 


Example 9 : 




lim 

x->+oo 



x 


Mathematica Code : 


F=Function[{x}, (^~) ]; 

Pl=Plot[{ArcTan[F[Tan[x]]]},{x,- n/2, n /2},PlotRange^>{0, n/2}]; 
P2=Plot[{ArcTan[e 2x2 {x,-n/2, n /2},PlotRange—>{0, n /2},PlotStyle^>Dashing[{.01}]]; 
P3=Plot[{ArcTan[F[Tan[x]]], ArcTan[e 2x2 ]}, {x,- n/2, n/2}, PlotRange—>{0, n/2}, PlotStyle n 
{Dashing[{l}], Dashing} {. 01}]}]; 

Show [ GraphicsArray[{Pl, P2, P3}]]; 
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As the Code and the Graphs show the result is e 2x2 . 
Checking the result : 

We know that : 

lim (1 + -) bx = e ab 

x-»+oo X 


So: 


li m x— >+oo 



x(l+f) 

x ^-f) 


X 


lim x ^ + 00 (l+|) x 

lim x _ +s 0 (l~) x 



e 2x - 2 


Example 10 : 

2x — 1 

lim 

x_> -°° V3x 2 + x + 1 


Mathematica Code : 


F=Function[{x},- 


2 x — 1 


■h 


/ V 3 x 2 +x+l " 

Pl=Plot[{ArcTan[F[Tan[x]]]},{x,-n/2, n /2},PlotRange^{- n/2, n/2}]; 


r —2. 


P2=Plot[{ArcTan[- F ]},{x,- n/2, n /2},PlotRange->{- n/2, n /2},PlotStyle->Dashing[{.01}]]; 

V3 


P3=Plot[{ArcTan[F[Tan[x]]],ArcTan[-p: ]},{x,- n/2, n /2},PiotRange -*{- n/2 , n/2}, 

V 3 

PlotStyle->{Thickness[.01 ],Dashing[{.01 }]}]; 



Show[ GraphicsArray[{Pl, P2, P3}]]; 



—2 

As the Code and the Graphs show the result is . 

V 3 

Checking the result: 

As we know : 


Vax n + bx n 1 + cx n 2 + ••• 
x -> ±oo 


~VS 



So: 

1 -™ 2x— 1 _ ] . 2x— 1 _ 2x— 1 _ ] 2x 2 

V3x2+x+1 - Um^.oo 

Example 11 : 

lim V x3 — 5x 2 + 2x + 1 — x 

X->oo 

Mathematica Code : 


F=Function[{x},Vx 3 — 5x 2 + 2x + 1 — x 7; 

Pl=Plot[{ArcTan[F[Tan[x]]]},{x,-n/2, n /2},P!otRange->{- n /2, n/2}]; 

P2=Plot[{ArcTan[-^-]},{x,- n/2, n /2},PlotRange->{- n/2, n /2},PlotStyle—>Dashing[{.01}]]; 
P3=Plot[{ArcTan[F[Tan[x]]],ArcTan[ ]},{x,-n/2, n/2},PlotRange -*{- n/2, n/2}, 
PlotStyle->{Thickness[.01 ],Dashing[{.01 }]}]; 

Show [ GraphicsArray[{Pl, P2, P3}]]; 



The graph and the code show that the result is -y . 

Checking the result : 

We know that : 

If n is a positive Integer number and p(x) is a polynomial in the following form 
p(x) = x n + biX 11-1 + b 2 x n “ 2 + ••• + b n _iX + b n 

Then 

ta([p(x)P-x)-| 


So 


3 1 —5 

lim -y x 3 — 5x 2 + 2x + 1 — x = — 

X^oo 3 


Comparison Between The Speed Of Different Functions 
In Growing To Infinity. 




/ t,s possible to compare the speed of different functions in growing to the Infinity using the 
"Infinity" plots . To do so we should plot different functions at the same time on the same 
graph and compare the speed of growing visually, 
here is some examples showing the method. 


Example : 

Comparison between the speed of (x,x 2 , e x , Ln(x)) in growing to infinity. 
Mathematica Code 


Fl=Function [{x},x]; 

F2=Function[{x},x A 2]; 

F3=Function[{x}, e A x]; 

F4=Function[{x},Log[x]]; 

Plot[{ArcTan[Fl[Tan[x]]],ArcTan[F2[Tan[x]]],ArcTan[F3[Tan[x]]] ,ArcTan[F4[Tan[x]]]},{x,- 

^/ 2 ,^/ 2 }]; 



As you can see , e x goes to +oo Faster than x 2 , and x 2 goes to+ °o faser than x . 

And Ln(x) goes to +oo slowly compared to (x,x 2 , e x ) 

Further Examples : 

We can use the above conclusion about the speed of functions to go to Infinity to determine the 
limits of functions. 


Example 1 : 

Try to find the the following limit : 


lim x — >+oo 


Ln(x)+x 2 

e x +x 


= ? 


Matematica Code : 


F=Function[{x},( Ln(x) + x A 2)/(e A x + x)]; 
Plot[ArcTan[ F[Tan[x]]], {x, - n/2,n/2 }] 



As you can see the function comes to zero as x -» +oo . 
Try evaluating the limit using Mathematica : 


Limit[- 


.Ln(x)+x 2 

e x +x 


, X -» +00] 


The result would be zero : 


0 


Example 2: 

Try finding the following limit : 

Mathematica Code : 


F=Function[{x}, (e A x) / (Ln(x)) ]; 
Plot[ArcTan[ F[Tan[x]]], {x, - n/2,n/2 }] 



Try evaluating the limit using Mathematica : 


+coJ 


The result would be oo; 


00 


As I mentioned before e x goes to +oo faster than Ln(x) . So the result is true . 


Asymptote 






nother helpful application of the Infinity Plots is to show the asymptote of functions 
visually . We can use the infinity plots to show the vertical , horizontal and inclined 
Asymptotes. We will show the method by some examples. 


Vertical Asymptote 

Vx — 2 

f 0) = — — 7 

x z — 4 


x 2 -4 = 0 


x = +2 


Butx = —2 is not Asymptote , because the Radical in the numerator of the function will have 
minus sign . 

Mathematica Code 


F=Fu notion [ 

Pl=Plot[{ArcTan[F[Tan[x]]]},{x,- n/2, n /2},PlotRange^>{- n /2, n/2}]; 

P2=Plot[{0},{x,0, n/2},GridLines^>{{ArcTan[2]},{0}},PlotRange^>{0, n/2}] 
P3=Plot[{ArcTan[F[Tan[x]]],},{x,0, n/2}, GridLines^>{{ArcTan[2]},{0}},PlotRange^>{0, n/2}, 
PlotStyle^>{Thickness[. 001], Dashing [{.01}]}]; 

Show] GraphicsArray[{Pl, P2, P3}]]; 



Horizontal Asymptote 

, s 5 — 3x 2 
f 0) = Y 

lim x ^ Too f(x) = 3 



Also we can check the Vertical Asymptote : 


, lim f(x) = — oo 

X->1 + 


lim f(x) = +oo 

X-+1 - 


, lim f(x) = —oo , lim f(x) = +oo 

X - > — 1 — x->— 1+ 


Math ematica Cock 

F= Fu notion [{x}, ^^7; 

Pl=Plot[{ArcTan[F[Tan[x]]]},{x,- n/2, n /2},PlotRange^>{- n /2, n/2}]; 
P2=Plot[{ArcTan[3]},{x,0, n/2},PlotRange^>{- n/2, n /2},PlotStyle->Dashing[{.01}]] 
P3=Plot[{ArcTan[F[Tan[x]]],ArcTan[3]},{x,- n/2, n /2},PlotRange^>{- n/2, n/2}, 
PlotStyle^>{Thickness[.001 ],Dashing[{.01 }]}]; 

Show[GraphicsArray[{Pl,P2,P3}]]; 





Inclined Asymptote 


f(x) = 


X J 


X 2 + 1 


y = ax + b a = lim 


f(x) 


x->oo x 


b = lim (f(x) — ax) 


a = lim 


f(x) 


= lim 


x->oo x 


x->oo x 2 + 1 


= 1 


b = lim (f(x) — ax) = lim 

X— >00 ” ' " 

-» y = x 


X J 


x->co X 2 + 1 


— x = lim 


—x 


x->co x 2 + 1 


= 0 


Mathematica Code 

„ X^ 

F=Function[{x},——]; 

X z + 1 

G=Func tion[{x},x]; 

Pl=Plot[{ArcTan[F[Tan[x]]]},{x,- n /2, n /2},PlotRange^>{- n/2, n/2}]; 

P2=Plot[{ArcTan[G[Tan[x]]]},{x,- n /2, n /2},PlotRange^>{- n/2, n/2}, PlotStyle — ► Dashing[{.01}]]; 
P3=Plot[{ArcTan[F[Tan[x]]],ArcTan[G[Tan[x]]]},{x,- n /2, n /2},GridLines->{2},PlotRange->{- n/2, 
n /2},PlotStyle->{Thickness[.001],Dashing[{.01}]}]; 

Show [ GraphicsArray[{Pl, P2, P3}]]; 

Limit [ F[x],x^>-lnfinity] 





Finding The Finite Number Of Real Roots Of Functions Visually 
Asxls Varying In The Rang — oo < x < +oo . 




/ 


/ the function does not hove an oscillating behavior as x ^ ±oo and if the number of roots 
of the function is finite , we can use the "Infinity" plot to visually determine the number of 
roots. Also it,s possible to use the Infinity transformation to find the roots . 


Example 1: 

We want to find the number of roots of the following function as x is varying in the range 

- oo < x < oo. 

sin(6 x (x 3 — .4)) + x 4 ) 

Mathematica Code 


F=Function[{x},Sin[6*(x A 3-.4)]+x A 4]; 

Plot[{ArcTan[F[Tan[x]]]},{x,-n/2,n/2}] 



From the plot we see that the function has 4 Real roots . 

Example 2: 

We want to find the number of roots of the following function as x is varying in the range 

- oo < x < oo. 


sin(12 * x) + x 4 

M a th ematica Code 

F= Function! {x},Sin[12 *x]+x A 4 ]; 
Plot[{ArcTan[F[Tan[x]]]},{x,-n/2,n/2}] 





From the plot we see that the function has 8 Real roots . 



53 


Finding The Finite Real Roots of Functions Visually using 
Infinity Plots. 




T lo visually find the finite real roots of the functions , first we should Infinity plot the 
function , then we should guess the roots using the Infinity Plot . Then we will use the 
Mathematica Find Root Command to find the roots through the transformed function. 
We will Introduce the method through an example . 


Example : 

Find the real roots of the following function : 


f = tan(2 x x) + x 3 
Solution : 

The Transformed function in Mathematica code will be : 


F=Function [x, Tan [2*x] +x A 3] ; 

G=Function [x, ArcTan [F [Tan [x] ] ] ] ; 

Plot [ArcTan[F [Tan [x] ] ] , { x, -Pi/2 , Pi/2 } ] 



From the picture we can guess the Real Roots approximated location . 


FindRoot [G [x] =0 , { x, { 0 , - . 81 , .81}}, Maxi ter at ions— >1000 ] 


Result : 



(x -» {0., -0.8350493275694328,0.835049327569433}} 



Real Roots : 


Tan[{0, -0.835049327569433,0.835049327569433}] 


Result : 


{0,-1.1045808419145513,1.1045808419145513} 
Testing the result using the main function : 


FindRoot[F[x] == 0, (x, {0, -1.1045808419145513,1.1045808419145513}}, Maxlterations 
-> 1000 ] 

Result : 


{x -> {0., -1.1045808419145513,1.1045808419145513}} 

So you can see that the result from the main function and from the transformed function using 
initial guesses optioned from the infinity plot is true . 
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• Developing A Code To Find The Finite Number Of Real 
Roots Of Functions As x Is Varying In The Range 

— 00 < x < +00 


II Tow we are going to develop a LUA Code to determine the finite number of Real roots of 
I II functions as x is varying in the range - oo < x < go . The procedure is that we will try to 
A. W find the finite number of Real roots of the transformation function in the range 
- j < x < j . Just note that this procedure is not suitable for alternative functions in other words 
we can not use functions with infinite number of roots in this method. 

Introduction to Lua Language : 

Lua is an extension programming language designed to support general procedural programming with 
data description facilities. It also offers good support for object-oriented programming, functional 
programming, and data-driven programming. Lua is intended to be used as a powerful, light-weight 
scripting language for any program that needs one. Lua is implemented as a library, written in clean C 
(that is, in the common subset of ANSI C and C++). 

Being an extension language, Lua has no notion of a "main" program: it only works embedded 
in a host client, called the embedding program or simply the host. This host program can invoke 
functions to execute a piece of Lua code, can write and read Lua variables, and can register 
C functions to be called by Lua code. Through the use of C functions, Lua can be augmented to 
cope with a wide range of different domains, thus creating customized programming languages 
sharing a syntactical framework. The Lua distribution includes a sample host program called lua, 
which uses the Lua library to offer a complete, stand-alone Lua interpreter. 

Lua is free software, and is provided as usual with no guarantees, as stated in its license. The 
implementation described in this manual is available at Lua's official web site, www.lua.org. 

Lua Codes: 

Here we are going to use the Kronecker - Picard theory to find the roots of functions . 

We have used the LNA "Lua Numerical Analysis " Package . 

What is LNA ? 

LNA stands for "Lua Numerical Analysis " ; it is a package of Numerical Analysis functions , 
together with several utility and plotting functions. All functions in the package are meant to be 
called from a program written in Cplua , which runs on a Casio Classpad 300 Calculator. 

Example : 

Find the roots of the following function as x ^ ±oo 
f(x) = sin(6x) + x 3 



Driver Program : 



require("LNAplot/PlotFunc","LNAplot/PlotData","LNA/Krone","LNAutils/Pi") 

local function g(x) 
return 

math.sin(6*x)+x A 3 

end 

local function f(x) 

return math.atan(g(math.tan(x))) 

end 

local function dfdx(x) return 

(6*math.cos(6*math.tan(x))*(l/math.cos(x)) A 2+3*(l/math.cos(x)) A 2*math.tan(x) A 2)/(l+(math. 

sin(6*math.tan(x))+math.tan(x) A 3) A 2) 

end 

local function d2fdx2(x) return 
(- 

18*(l/math.cos(x)) A 4*(2*math.cos(6*math.tan(x))+math.tan(x) A 2) A 2*(math.sin(6*math.tan(x)) 
+math.tan(x) A 3)+6*(l/math.cos(x)) A 2*(2*math.cos(6*math.tan(x))*math.tan(x)+math.tan(x) A 3 
+( 1/math. cos(x)) A 2*(- 

6*math.sin(6*math.tan(x))+math.tan(x)))*(l+(math.sin(6*math.tan(x))+math.tan(x) A 3) A 2))/(l+( 

math.sin(6*math.tan(x))+math.tan(x) A 3) A 2) A 2 

end 

local Nr, roots, axesdata 

roots, Nr=KroneRoots(f,dfdx,d2fdx2,-Pi/2,Pi/2) 
print(Nr.." roots found:") 

for i,v in ipairs(roots) do print(i,math.tan(v)) end 
froots={) 

for j=l,Nr do froots[j]=g(math.tan(roots[j])) end 
print("\nFunction values at roots:") 
for i,v in ipairs(froots) do print(i,v) end 

axesdata=PlotFunc(f,{-Pi/2,Pi/2),{-Pi/2, Pi/2 }, false) 

PlotData(roots,froots,{),axesdata,true,2) 


Screen Capture Showing The Results : 


V Break 


5 roots found: 

1 -0. 90689035948226 

2 -0. 55172257558493 

3 0 

4 0.55172257553493 

5 0. 90609035940226 

Function values at roots: 

1 1.1 102230246252e-016 

2 -1 . 30777070070 1 4e-0 1 6 

3 0 

4 1 . 30777070070 14e-0 16 

5 -1 . 1 102230246252e-016 
Done. 

Press [EXE]... 


Done. frill 


Screen Capture Showing The Results Visually: 



Checking The Results Using Mathematica : 

Mathematica Code 


FindRoot[Sin[6*x]+(x) A 3==0 ' {x, {-.9, -.5,0, .5,. 9}}] 


Result : 


{x— >{-0.90689, -0.551723,0. ,0.551723,0.90689}} 



As you can see the results are true. 


Requirements Pi : 

Please note that the following documents ( the requirement section ) are Copyrighted by PAP 
[Developer oF the LNA (Lua Numerical Analysis) Package] . 

KroneRoots 

The function KroneRoots computes all the roots of a function within a given interval , or , 
optionally , a prescribed number of roots . It can be used locating and computing all roots of a 
function within a given interval .The algorithm implemented in this function is based on the 
Kronecher-Picard theory (hence its name) , and uses recursive auxiliary functions. It is currently 
the most complex algorithm included in LNA .Despite its complexity , however ; the function can 
be easily called in a user program. 

Syntax 


roots ,Nr=KroneRoots(f,dfdx,d2fdx2,a,b,r_reg,xi) 


returns a vector roots , containing all the roots of the function f inside the interval {. a,b }; 
optionally , it may return only a prescribed number of roots . This function also returns the total 
number of roots , Nr. The arguments dfdx, and d2fdx2 define the first and second derivatives 
of the function , respectively. The arguments r_reg controls how many roots should be 
returned. If omitted , or if r_reg<= 0 , all the roots will be returned . The optional argument xi 
defines an appropriate number, such that xi*dfdx(x) is not too small ( or too large) compared 
to f(x),for all x inside the interval {a,b} . The default value is xi=l. 

Bisect 

The function Bisect computes the root of a function within a given interval via the Bisection 
method .This is the simplest ( and the slowest) method for root finding . However, its 
convergence is always guaranteed , and it is ofen used by more complex numerical methods . 

The algorithm implemented in Bisect is a modification of the "classic" algorithm , sometimes 
called "Simplified Bisection", and it is slightly faster. 

Syntax 


root, error=Bisect(i f,xl,xr, eps,maxit,sho w) 


returns the root of the function f inside the interval {. xl,xr} , together with estimation of the 
absolute error. The arguments eps, maxitand show are optional. Eps defines the desired 
accuracy (default : Epsilon, i.e., 1.12xl0' 16 ); it can be set to "auto", which is equivalent to the 
default accuracy, maxit defines the maximum number of iterations ( default : 100) . show is a 


Boolean argument that control wether progress of the iterative process will displayed or not 
(default : false); if set to true, the function value , f(root) , at each iteration will be displayed. 


Romberg : 

The function Romberg computes the definite integral of a function via the Romberg Method. 
Syntax : 

q=Romberg(f,a,b,eps,k,show) 


returns the definite integral ,q, of the function f, integrated from ato b. The arguments eps 
,/r,and show are optional : epsdefines the desired accuracy . ^defines the order of the method 
; show\s a Boolean argument that controls wether an integration progress will be displayed or 
not. 


LNA Utility Constants And Functions : 

Epsilon 

The constant Epsilon defines the smallest positive number, e , which satisfies inequality 1+ £ > 1 .Due 
to computer arithmetic , adding a very small number to unity result exactly one. This constant defines 
Epsilon as s = 1.12810" 16 


Pi 

The constant Pi defines ti . Since ji is often used in numerical computations , it is defined as a 
LNA constant , accurate to 16 digits . This constant defines Pi as ti=3. 1415927410125732 


Part 

The function Part return part of a vector A, containing elements from Afimin] to A[imax] 

Syntax 

P=Part(A,imin.imax) 


Return a vector P , with imax-imin+1 elements , so that P[l]=A[imin] , P[2]=A[imin+l] , and so 
on, until P[imax-imin+l]=A[imax] 


The codes required to run the driver program PI 



The Following codes are copyrighted by PAP [Developer of LNA ( Lua Numerical Analysis ) 
Package] 


LNA/Krone 


requiref "table", "LN A/Bisect", "LNA/Romberg") 

local f clfdx, d2fdx2 
local roots 
local xi=l 
local nr,cr 

local function Kronelntg(x) 
local dfsq=dfdx(x) A 2 

return (f(x) *d2fdx2(x)-dfsq)/(f(x) A 2+xi A 2 *dfsq) 
end 

local function Nroots(a,b) 
local Nr 

Nr=Romberg(Kronelntg, a, b,0.1) 

Nr=-0.3183098861837907*(xi*Nr-math.atan(xi*dfdx(b)/f(b))+math.atan(xi*dfdx(a)/f(a))) 

Nr=math.floor(Nr+0.49999999999999) 

return Nr 

end 

local function Roots(a,b,Nr) 
local subs,h,xl,xr,Njr 
if Nr==l then 
cr=cr+l 

roots[cr]=Bisect(f a,b) 
else 

subs=math.max(2,Nr) 
h=(b-a)/subs;xr=a 
for j=l,subs do 
xl,xr=xr,xr+h 
if j<subs then 
Njr=Nroots(xl,xr) 
else 

Njr=Nr-Njr 

end 

if Njr>0 then 



Roots(xl,xr,Njr) 
if cr==nr then 
return 
end 
end 
end 
end 
end 

function KroneRoots(g, dgdx, d2gdx2, a,b,...) 
f=g;dfdx=dgdx;d2fdx2=d2gdx2 
local Nr 
if a==b then 

print("Error: KroneRoots: Wrong inteval") 
elseif a>b then 
a,b=b,a 
end 

if a rg["n"]>=2 then 
xi=arg[2] 
end 

Nr=Nroots(a,b) 
if Nr>0 then 
if arg["n"]>=l then 
ifarg[l]<=0 then 
nr=Nr 

elseif arg[l]<=Nr then 
nr=arg[l] 
else 
nr=Nr 

print("Warning: KroneRoots: Interval contains fewer roots than those required.") 
end 
else 
nr=Nr 
end 

roots={} 

cr=0 

Roots(a,b,Nr) 

else 

print("Warning: KroneRoots: No simple roots found.") 
end 

return roots, Nr 
end 

export{KroneRoots=KroneRoots) 



LNA/Bisect 


requiref "LNA utils/E psilon ") 

function Bisect(fxl,xr,...) 
local eps=Epsilon 
local maxit=100 
local show=false 
if arg["n"]>=l then 
if arg[l ]~="auto " then 
eps=arg[l] end 
if arg["n "]>=2 then 
maxit=arg[2] 
if arg["n"]>=3 then 
show=arg[3] 
end 
end 
end 

local root, error 
local xp,xc,fc,sl,sp,half 
xp=xl;xc=xl 
sl=math.sign(l,f(xl)) 
sp=sl;half=xr-xl 
if sl*math.sign(l,f(xr))>0 then 

print(" Error: Bisect: Interval does not contain a root, or contains an even number of roots. ") 
return 
end 

if show then print(" i discrepancy") end 

for i=l,maxit do 

half=0.5*half 

xc=xp+sl*sp*half 

fc=f(xc) 

if show then printf("%2i %+e\n",i,fc) end 
if half <=eps then 
root=xc;error=h alf 
return root, error 
end 

iffc==0 then 
root=xc;error=0 
return root, error 
else 
xp=xc 

sp=math.sign(l,fc) 

end 
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end 

root=xc;error=half 

print("Warning: Bisect: Iterations limit reached ' but the accuracy is not satisfied. ") 

return roof error 

end 

export{Bisect=Bisect} 


LNAutils/Epsilon 


Epsilon=1.12E-16 

export{Epsilon=Epsilon} 


LNA/Romberg 


requireftable ", "LNAutils/Part") 

local function TrapStepff a,b, t,n) 
local s 

local it,dei,x,su 
if n==l then 
s=0.5*(b-a)*(f(a)+f(b)) 
else 

it=2 A (n-2);del=(b-a)/it 

x=a+0.5*del 

su=0 

for j=l,it do 
su=su+f(x);x=x+del 
end 

s=0.5*(t+(b-a)*su/it) 

end 

return s 
end 

local function Pollntrp(xa,ya,x) 
local y,dy 

local n=table.getn(xa) 
local ns,difdift,c,d,ho 
,hp,w,den 
c=table.copy(ya) 
d=table.copy(ya) 



ns=l 

dif=math.abs(x-xa[l ]) 
for i=l,n do 
dift=math . abs(x-xa[i ] ) 
if dift<dif then 
ns=i;dif=dift 
end 
end 

y=ya[ns];ns=ns-l 
for m=l,n-l do 
for i=l,n-m do 
ho=xa[i]-x 
hp=xa[i+m]-x 
w=c[i+l]-d[i] 
den=ho-hp 
if den ==0 then 

print(" Error: Romberg: Roundoff error. ") 
return 
end 

den=w/den 

d[i]=hp*den 

c[i]=ho*den 

end 

if 2* ns<n-m then 
dy=c[ns+l] 
else 

dy=d[ns];ns=ns-l 

end 

y=y+dy 

end 

return y,dy 
end 

local function Reached(q,dq,eps) 
local ok 
if eps>0 then 
ok=math.abs(dq)<eps 
else 

ok=math.abs(dq)< 

-eps*math.abs(q) 

end 

return ok 
end 

function Romberg(fa,b,...) 
local q 


local eps=lE-6 
local k=2 
local show= false 
local km,dq,h,s 
local jmax=20 
local jmaxp=jmax+l 
if org["n"]>=l then 
eps=arg[l] 
if arg["n"]>=2 then 
k=arg[2] 

if arg["n"]>=3 then 
show=arg[3] 
end 
end 
end 
km=k-l 
h={};s={} 
h[l]=l 

for j=l,jmax do 
s[j]=TrapStep(f a , b,s[j],j) 
if j>=k then 

q,dq=Pollntrp(Part(h,j-km,j),Part(s,j-km,j),0) 

if show then print(j,q) 
end 

if Reached(q,dq,eps) then 

return q 
end 
end 

s[j+l]=s[j] 

h[j+l]=0.25*h[j] 

end 

print("Warning: Romberg: Maximum number of iterations reached, but the acuracy criterion is 
not satisfied. ") 
return q 
end 

export{Romberg=Romberg} 


LNAutils/Part 


function Part(A,imin,imax) 

local P 

P={} 


/ 


for i=imin,imax do 
P[i-imin+l]=A[i] 
end 

return P 
end 

export{Part=Part} 


PlotFunc 


require ("table", "draw", "LNAplot/PlotUtil") 

function PlotFunc(f,xv,yv,...) 

-Optional arguments: 
local wait=true 
local lwidth={l} 
local tics={" auto", "auto"} 
local grid=true 
local lblpos={0,0} 
local lblsize={9,9} 
local c=" auto" 
local discont={} 

local scale 

local lws,funcs,x,xp,X,Xp, Y, Yp,disconts 
local id={ } 
if arg["n"]>=l then 
wait=arg[l] 
if arg["n"]>=2 then 
lwidth=arg[2] 

if type(lwidth)=="number" then 
lwidth={lwidth} 
end 

if a rg["n"]>=3 then 
tics=arg[3] 

iftype(tics)~="table" then 
tics={ticc,tics} 

end 

ifarg["n"]>=4 then 
grid=arg[4] 

if arg["n "]>=5 then 
lblpos=arg[5] 

iftype(lblpos)=="number" then 


lblpos={ Iblpos, Iblpos} 

end 

ifarg["n"]>=6 then 
lblsize=arg[6] 

if type(lblsize)=="number" then 
lblsize={lblsize,lblsize} 

end 

if arg[”n"]>=7 then 
c=arg[7] 

if arg["n"]>=8 then 
discont=arg[8] 

end 

end 

end 

end 

end 

end 

end 

end 

iws=#iwidth 

if type(f)== "function " then 

f={f} 

end 

funcs=#f 

for i=lws+l,funcs do 
lwidth[i]=lwidth[lws] 
end 

discont,disconts=Discontinuities(funcs,discont) 
for i=l,funcs do 

id[i]=l 

table.sort(discont[i]) 

end 

draw.onbuffer() 
ifxv[l]~=nil then 
scale=SetScale(xv,yv) 

PlotAxesfxv, yv, scale, tics, grid, Iblpos, I b I size, c) 
else 

xv=yv[l] 

scale=yv[3] 

yv=yv[2] 

end 

Xp=0;xp=xv[l] 

Y=0;Yp={} 
for i=l,funcs do 
Yp[i]=(yv[2]-f[i](xv[l]))*scale[2] 
end 


for X=l,X_pixels do 
x=xv[l ]+X/scale[l ] 
for i=l,funcs do 
Y[i]=(yv[2]-f[i](x))*scale[2] 

if id[i]>disconts[i] or x<discont[i][id[i]] or xp>discont[i][id[i]] then 
draw. linefXp, Yp[i],X , Y[i], 1, /width [i]) 
elseifxp<=discont[i][id[i]] and x>=discont[i][id[i]] then 
id[i]=id[i]+i 
end 

Yp[i]=Y[i] 

end 

xp=x;Xp=X 

end 

showgraphf) 
draw.update() 
if wait then waitkeyf) end 
showconsolef) 
return {xv,yv, scale} 
end 

export{PlotFunc=PlotFunc} 


LNAplot/PlotUtil 


require ("string", "table", "draw", " LNAutils/OrderMag ") 

X_max=158 

Y_max=213 

X_center=79 

Y_center=106 

X_pixels=159 

Y_pixels=214 

local function SetScale(xv,yv) 
local xscale,yscale 
xscale=X_pixels/(xv[2]-xv[l ]) 
yscale=Y_pixels/(yv[2]-yv[l]) 
return {xscale,yscale} 
end 

local function AutoXtics(range) 
-tics=math.max(math.floor( range/7), 1) 
local tics=OrderMag(range) 



if range<4*tics then 
tics=tics/2 
end 

return tics 
end 

local function AutoYtics(range) 
~tics=math.max(math.floor((yv[2]-yv[l])/9),l) 
local tics=OrderMag(range) 
if range<4*tics then 
tics=tics/2 
end 

return tics 
end 

function Discontinuities(funcs,discont) 
local D=table.copy(discont) 
local id={ } 
local disconts={} 
for i=l,funcs do 
if D[i]==nil then 
D[i]={ } 

el seif type(D[i])=="number" then 
D[i]={D[i]} 
end 

discon ts[ i }=#D[ i ] 
end 

return D,disconts 
end 

function PlotAxes(xv,yv,scale,tics,grid,lblpos,lblsize,c) 
local Xc, Yc,X, Y,XI, Yl,lblxp,lblyp 
local lblcut=lE-8 
-Auto tics selection: 
if tics[l]=="auto" then 
tics[l ]=AutoXtics(xv[2]-xv[l]) 
end 

if tics[2]=="auto" then 
tics[2]=AutoYtics(yv[2]-yv[l]) 
end 

if c==" auto" then 
c={0, 0} 
end 

—[[ AXES ]]-- 
Xc=(c[l]-xv[l ]) *scale[l ] 



Yc=(yv[2]-c[2])*scale[2] 
draw. line(0, Yc,X_max, Yc) 
draw. line(Xc,0,Xc, Y_max) 

-[[ TICS & GRID //- 
if lblpos[l]==0 then 
YI=Y_pixels-lblsize[l ] 
else 
YI=Yc+l 
end 

if lblpos[2]==0 then 
Xl=2 
else 

XI=Xc+2 

end 

for xp=c[l ]-math.floor((c[l]-xv[l ])/tics[l ])*tics[l ],math.floor((xv[2]-c[l ])/tics[ 1]) *tics[ 1 ], tics[ 1 ] 
do 

X=(xp-xv[ 1 ])*scale[ 1 ] 
draw. line(X, Yc-2,X, Yc+2) 
if grid then 
for Y=0, Y_pixels,4 do 
draw.pixel(X,Y) 
end 
end 

if I blsize [1]>0 and X>0 and X<X_pixels then 
lblxp=xp 

if math.abs(xp)<=lblcut then 
lblxp=0 
end 

draw. text(X+2, Yl,lblxp,l,lblsize[l ]) 
end 
end 

for yp=c[2]-math.floor((c[2]-yv[l])/tics[2])*tics[2],math.floor((yv[2]-c[2])/tics[2])*tics[2],tics[2] 
do 

Y=(yv[2]-yp) *scale[2] 
draw.iine(Xc-2, Y,Xc+2, Y) 
if grid then 
for X=0,X_pixeis,4 do 
draw.pixei(X,Y) 
end 
end 

if I blsize [2] >0 and Y>0 and Y<Y_pixels then 
lblyp=yp 

if math.abs(yp)<=lblcut then 
lblyp=0 
end 

draw. text(XI, Y-lblsize[2],lblyp / l / lblsize[2]) 



end 

end 

end 

function PlotPointfX, Y,pointtype, pointsize) 

-Non-filled circle: 
if pointtype==0 then 
draw.point(X, Y,l, pointsize) 
elseif pointtype==l then 
dra w. circle(X, Y, pointsize, 1,1,-1) 
draw.point(X,Y) 

-Crossed circle: 

elseif pointtype==2 then 

dra w. cirde(X, Y, pointsize , 1,1,-1) 

draw. line(X-pointsize, Y,X+ pointsize, Y) 

draw. line(X, Y-pointsize,X, Y+pointsize) 

-Filled circle: 

elseif pointtype==3 then 

draw. cirde(X, Y, pointsize, 1,1, 1) 

-Non-filled rectangle: 
elseif pointtype==4 then 

draw. rect(X-pointsize, Y-pointsize,X+pointsize, Y+pointsize, 1, 1, -1 ) 
draw.point(X,Y) 

-Crossed rectangle: 
elseif pointtype==5 then 

draw. rect(X-pointsize, Y-pointsize,X+pointsize, Y+pointsize, 1, 1, -1 ) 
draw. line(X-pointsize, Y,X+pointsize, Y) 
draw.linefX, Y-pointsize,X, Y+pointsize) 

-Filled rectangle: 
elseif pointtype==6 then 

draw. rect(X-pointsize,Y-pointsize,X+pointsize, Y+pointsize, 1,1,1) 

end 

end 

local function DrawLabelf) 
return 
end 

export{X_max=X_max, Y_max=Y_max,X_cen ter=X_ center, Y_center=Y_center,X_pixels=X_pixels, Y 
_pixels=Y_pixels,SetScale=SetScale,Discontinuities=Discontinuities,PlotAxes=PlotAxes,PlotPoint=P 
lotPoint} 


LNAutils/OrderMag 
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require("string") 
function OrderMag(x) 

local xexp=string.lower(string.format("%e",x)) 
local e=string.find(xexp, "e") 
local n=string.len(xexp) 
return 10 A string.sub(xexp, e+l,n) 
end 

export{ OrderMag=OrderMag } 


LNAplot/PlotData 


require ("table", "draw", "LNAplot/PlotUtil") 

function PlotData(x, y,xv,yv,...) 

-Optional arguments: 

local wait=true 

local ptype={l} 

local psize={3} 

local lwidth={0} 

local tics={" auto", "auto"} 

local grid=true 

local lblpos={0,0} 

local lblsize={9,9} 

local c="auto" 

local points=#x 
local scale 

local sets,pws,lws,X,Xp,Y,Yp 
if arg["n"]>=l then 
wait=arg[l] 
if arg["n"]>=2 then 
ptype=arg[2] 

if type(ptype)==" number" then 
ptype={ptypej 
end 

if a rg["n"]>=3 then 
psize=arg[3] 

if type(psize)==" number" then 
psize={psize} 

end 

if arg["n"]>=4 then 



Iwidth=arg[4] 

if type(lwidth)=="number" then 
lwidth={lwidth} 

end 

if arg["n"]>=5 then 
tics=arg[5] 

iftype(tics)~="table" then 
tics={tics,tics} 

end 

ifarg["n"]>=6 then 
grid=arg[6] 

if org["n "]>=7 then 
lblpos=arg[7] 

if type(lblpos)=="number" then 
lblpos={lblpos , Iblpos } 

end 

ifarg["n"]>=8 then 
lblsize=arg[8] 

if type(lblsize)=="number" then 
lblsize={lblsize,lblsize} 

end 

if arg["n"]>=9 then 

c=arg[9] 

end 

end 

end 

end 

end 

end 

end 

end 

end 

if type(y[l])=="number" then y={y} end 
sets=#y 
pws=#ptype 
for i=pws+l,sets do 
ptype[i]=ptype[pws] 
end 

pws=#psize 
for i=pws+l,sets do 
psize[i ]=psize[ pws] 
end 

lws=#lwidth 

for i=lws+l,sets do 

lwidth[i]=lwidth[lws] 

end 



draw.onbuffer() 
ifxv[l]~=nil then 
scale=SetScale(xv,yv) 

PlotAxes(xv,yv,scale, tics, grid, Iblpos, lblsize,c) 
else 

xv=yv[l] 

scale=yv[3] 

yv=yv[2] 

end 

Xp=(x[l ]-xv[l])*scale[l] 

Y=0;Yp={} 
for j=l,sets do 
Yp[j]=(yv[2]-y[j][l]) *scale[2] 
PlotPoint(Xp,Yp[j],ptype[j],psize[j]) 
end 

for i=2,points do 
X=(x[i]-xv[ 1]) *scale[ 1 ] 
for j=l,sets do 
Y[j]=(yv[2]-y[j][i]) *scale[2] 
if psize[j]>0 then 
PlotPoint(X,Y[j],ptype[j],psize[j]) 
end 

if lwidth[j]>0 then 
draw. line(Xp, Yp[j],X, Y[j],l,lwidth[j]) 
Yp[j]=Y[j] 
end 
end 
Xp=X 
end 

showgraph() 
draw.update() 
if wait then waitkeyf) end 
showconsolef) 
return {xv,yv, scale} 
end 

export{PlotData=PlotData} 


LNAutils/Pi 


Pi=3. 1415927410125732 



export{Pi=Pi} 



Developing A Method To Find The Finite Number Of Infinities 
Of The Functions . 




7 o visually find the finite number of Infinities of the functions, first we should Infinity 
plot the function , then we should guess the Infinities using the Infinity Plot .Then we 
will use the Mathematica Find Root Command to find the Infinities through the 
Transformed function . We will Introduce the method through examples . 


Example 1: 

Find the Infinites of the following function : 


5 



Solution : 

The Transformed function in Mathematica code will be : 


ArcTan[- 


5 + Tan[x] 


r] 


The Infinity plot will be : 


Plot[ArcTan[- 


5 n n 

-5 + Tan[x] 2 ^ X, ~2 , 2^ 




Using the plot we can guess the Infinities . +oo (Positive Infinities) of the function can be find with the 
Initial guesses of {—1.16,1.16} through the transformed function . 

The Mathematica code to find the Positive Infinities of the transformed function is as follows : 


FindRoot[-| + ArcTan[- ? - p l^ M2 


],0, {-1.16,1.16}}] 


Result : 


{x -> {-1.1502619915113772', 1.1502619915109316'}} 


The Positive Infinities of the main function will be find using the following method (tangent of the results 

): 


Tan[{-1.1502619915113772', 1.1502619915109316'}] 


{-2.236067977502464,2.2360679774997902} 



Result : 



—oo (Negative Infinities) of the function can be find with the Initial guesses of (—1.15,1.15) through the 
transformed function . 

The Mathematica code to find the Negative Infinities of the transformed function is as follows : 


n 


FindRoot[+ — + ArcTanf- 


—5 + Tan[x] 


r ], (x, {-1.15,1.15}}, AccuracyGoal -> 7] 


Result : 


(x -> (-1.1502619845593214", 1.1502619753619217"}} 


The Negative Infinities of the main function will be find using the following method (tangent of the 
results ): 


Tan[{— 1.1502619845593214", 1.1502619753619217"}] 


Result : 


{-2.2360679357901296,2.2360678806057344} 


Testing the results : 

To test the results we will use the general methods of finding the asymptotes of the functions . The 
vertical Asymptote of the function will be find through the roots of the denominator of the function . 

Vertical Asymptotes : 

The Mathematica code to find the vertical Asymptotes is as follows : 

A[Solve[x 2 — 5 — 0,x]] 

Result : 

{{x — > -2.23606797749979"}, (x-> 2.23606797749979"}} 

You can see that the results are true . 



Example 2 : 


Find the Infinites of the following function : 


5 

^ (x — 1) x (x + 3) x (x — 6) 

Solution : 

The Transformed function in Mathematica code will be : 


ArcTan[(-_6 + Tan[x])(— 1 + Tan[x])(3 + Tan[x])^ 


The Infinity plot will be : 


Plot 


Arc T *n[ ( . 6 + t^,, + TanM)(3 + ^ ~ ^ 


The Result will be as follows : 



Using the plot we can guess the Infinities . +oo (Positive Infinities) of the function can be find with the 
Initial guesses of {—1.2, 0.75, 1.42} through the transformed function . 




The Mathematica code to find the Positive Infinities of the transformed function is as follows : 



n 


FindRoot[-- 


+ ArcTan[ ( ,_ 6 + Tan ^^_ 1 + Tan [ x ])( 3 + Tan[x]) 
-» 150, AccuracyGoal -» 6] 


], { x , (—1.2,0.75,1.42)}, Maxlterations 


Result : 


(x -> (-1.2490457631598328", 0.7853981415952962', 1.40564765150521791} 


The Positive Infinities of the main function will be find using the following method (tangent of the results 

): 


Tan[{-1.249045763 1598328", 0.7853981415952962", 1.4056476515052179'}] 


Result : 


{— 3.,1.,6.} 


or the result with 15 digits after the point will be : 


(-2.999999907615787", 0.9999999563956967", 6.0000000786230805"} 


—oo (Negative Infinities) of the function can be find with the Initial guesses of(— 1.3,0.81,1.39} through 
the transformed function . 


The Mathematica code to find the Negative Infinities of the transformed function is as follows : 


FindRoot[+ — 

5 

+ ArcTan[^_6 + Tan[x])(— 1 + Tan[x])(3 + Tan[x]) 
— * 4] 


],{x, (-1.3,0.81,1.39}}, 


AccuracyGoal 


Result : 


(x -> (-1.2490468690695897", 0.7854009174907848", 1.4056475158794814"}} 



The Negative Infinities of the main function will be find using the following method (tangent of the 
results ): 


Result : 


Tan[{-1.2490468690695897', 0.7854009174907848", 1.4056475158794814 s }] 


{-3.00001,1.00001,6.} 

or the result with 15 digits after the point will be : 

{-3.000010966749434,1.000005508201843,5.999995060474788) 

Testing the results : 

To test the results we will use the general methods of finding the asymptotes of the functions . The 
vertical Asymptote of the function will be find through the roots of the denominator of the function . 

Vertical Asymptotes : 

The Mathematica code to find the vertical Asymptotes is as follows : 

N[Solve[(x — 1) * (x + 3) * (x — 6) == 0,x]] 

Result : 

{{x -> -3. }, {x -> 1. }, {x -» 6. }} 

You can see that the results are true . 

Example 3 : 

Find the Infinites of the following function : 

5 

^ Sin[6 * (x 3 — 0.4)] + x 4 
Solution : 

The Transformed function in Mathematica code will be : 

5 

ArcTan[ (Tan[x] 3 _ p 4)] _|_ Tan[x] 4 ^ 


The Infinity plot will be : 


Plot[ArcTan{ 


Sin[6(— 0.4 + Tan[x] 3 )] + Tan[x] 


n n n n 

M x >~n’ nl PlotRange -»{--, -}] 


2 2 


2 2 


The Result will be as follows : 



Using the plot we can guess the Infinities . +oo (Positive Infinities) of the function can be find with the 
Initial guesses of {—0.8, —0.73, —0.48,0.63} through the transformed function . 

The Mathematica code to find the Positive Infinities of the transformed function is as follows : 


n 

FindRoot[— — 

5 

+ A rc Tan[^. n |.^_Q ^ + j an |y|3)] _|_ Tan[x] 
-» 150] 


-j], {x, {—0.8, —0.73, —0.48,0.63}}, Maxlterations 


Result : 


{x -» {-0.7789523837360715, -0.7442419437453658, -0.4518738358393392, 
0.6171284554985426}} 


The Positive Infinities of the main function will be find using the following method (tangent of the results 

): 


Tan[(-0.7789523837360715, -0.7442419437453658, -0.4518738358393392, 
0.6171284554985426}] 


Result : 


(-0.9871908283758791,-0.9208984161497941,-0.48536824294300907, 
0.7095827903029706} 


—oo (Negative Infinities) of the function can be find with the Initial guesses of 
(—0.77, —0.75, —0.44,0.6} through the transformed function . 

The Mathematics code to find the Negative Infinities of the transformed function is as follows : 
n 

FindRoot[+ — 

+ ArcTa„[ sin[6( _ 04 + T J M 3 )]+TanM J. {*. {-0.77, -0.75, -0.44,0.6}}, Maxlteratlons 
-» 150, AccuracyGoal -> 6] 


Result : 


{x -> (-0.7789522798492675, -0.7442420311869101, -0.45187363699060945, 
0.6171281837619244}} 


The Negative Infinities of the main function will be find using the following method (tangent of the 
results ): 


Tan[(-0.7789522798492675, -0.7442420311869101, -0.45187363699060945, 
0.6171281837619244}] 


Result : 


(-0.9871906232466549, -0.9208985777464935, -0.4853679972490557, 
0.7095823817449417} 


Testing the results : 

To test the results we will use the general methods of finding the asymptotes of the functions . The 
vertical Asymptote of the function will be find through the roots of the denominator of the function . 

Vertical Asymptotes : 

The Mathematica code to find the vertical Asymptotes is as follows : 


FindRoot[Sin[6 * (x 3 — 0.4)] 

+ x 4 , {x, {-0.9871908283758791, -0.9208984161497941, -0.48536824294300907, 
0.7095827903029706})] 


Result : 


{x -> {-0.9871908283758782, -0.9208984161497942, -0.48536824193310446, 
0.7095827903029706}) 


You can see that the results are true . 

Lua Driver Program to find the Vertical Asymptotes are as follows : 

We should break the Domain to 3 domains . The Domains are : 

n n 

— — < X < 1 , — 1 < X < 1 , 1 < x < — 


require("LNAplot/PlotFunc", "LNAplot/PlotData", "LN A/Krone", "LNAutils/Pi") 

local function g(x) 
return 

math.sin(6*(x A 3-0.4))+x A 4 

end 

local function f(x) 
return 

math.atan(math.sin(6*(-0.4+(math.tan(x)) A 3))+(math.tan(x)) A 4) 

end 

local function dfdx(x) return 
(18*math.cos(6*(- 

0.4+(math.tan(x)) A 3))*(l/math.cos(x)) A 2*(math.tan(x)) A 2+4*(l/math.cos(x)) A 2*(math.tan(x)) A 3 

)/(l+(math.sin(6*(-0.4+(math.tan(x)) A 3))+(math.tan(x)) A 4) A 2) 

end 

local function d2fdx2(x) 
return 

-((2 *(18*math.cos( 6 *( - 

0.4+(math.tan(x)) A 3))*(l/math.cos(x)) A 2*(math.tan(x)) A 2+4*(l/math.cos(x)) A 2*(math.tan(x)) A 3 

) A 2*(math.sin(6*(-0.4+(math.tan(x)) A 3))+(math.tan(x)) A 4))/(l+(math.sin(6*(- 

0.4+(math.tan(x)) A 3))+(math.tan(x)) A 4) A 2) A 2)+(36*math.cos(6*(- 

0.4+(math.tan(x)) A 3))*(l/math.cos(x)) A 4*math.tan(x)+12*(l/math.cos(x)) A 4*(math.tan(x)) A 2+3 

6*math.cos(6*(- 

0.4+(math.tan(x)) A 3))*(l/math.cos(x)) A 2*(math.tan(x)) A 3+8*(l/math.cos(x)) A 2*(math.tan(x)) A 4 


-324*(l/math.cos(x)) A 4*math.sin(6*(-0.4+(math.tan(x)) A 3))*(math.tan(x)) A 4)/(l+(math.sin(6*(- 

0.4+(math.tan(x)) A 3))+(math.tan(x)) A 4) A 2) 

end 

local Nr, roots, axesdata 

roots, Nr=KroneRoots(f, dfdx, d2fdx2, -1,1) 
print(Nr.." roots found:") 
if Nr~=0 then 

for i,v in ipairs(roots) do print(i,math.tan(v)) end 
froots={} 

for j=l,Nr do froots[j]=g(math.tan(roots[j])) end 
print("\nFunction values at roots:") 
for i,v in ipairs(froots) do print(i,v) end 

axesdata=PlotFunc(f {-Pi/2, Pi/2 }, {-Pi/2, Pi/2 /false) 

PiotData(roots,froots,{), axesdata, true, 2) 
end 


The Output for the specified domains will be as follow : 

--< x < 1 — 1 < x < 1 1 < x < - 

2 2 


Console Break 


Console Break 


Console Break 

Warning: KroneRoots: No 
simple roots found. 

0 roots found in the Domain 
-Pi/2<x<-l: 

Done. 

Press [<-] or [EXE]... 


4 roots found: 

1 -0.98719082837588 

2 -0. 92089841814979 

3 -0.4853682419331 

4 0.70958279030297 

Function values at roots: 

1 -2. 442490654 1 753e-0 1 5 

2 -3. 330669073S755e-016 

3 5. 2735593669695e-016 

4 -6. 6613331477509e-016 
Done. 

Press [<-] or [EXE]... 


Warning: KroneRoots: No 
simple roots found. 

0 roots found in the Domain 
l<x<Pi/2: 

Done. 

Press [<-] or [EXE]... 

Done. ijm] 


Done. ijm] 

Done. ijm] 
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Example 4 : 

Find the Infinites of the following function : 
f = tan(2 x x) + x 3 
Solution : 

The Transformed function in Mathematica code will be : 
ArcTan[Tan[2 * 7an[x]] + Tan[x] 3 ] 


The Infinity plot will be : 


Plot[^rcTan[Tan[2 * Tan[x]] + Tan[x] 3 ],{x, — — },PlotRange All] 


The Result will be as follows : 





Using the plot we can guess the Infinities . +oo (Positive Infinities) of the function can be find with the 
Initial guesses of {-1.3216,-1.17, —0.7, 0.6,1) through the transformed function . 

The Mathematica code to find the Positive Infinities of the transformed function is as follows : 


FindRoot[— 7t/2 + ArcTan[Tan[2 * Tan[x]] 

+ Tan[x] A 3], { x , {—1.3216, —1.17, —0.7, 0.6, 1,1. 3}}, Maxlterations 
-> 1000, AccuracyGoal -> 6] 


Result : 


{x -> {-1.321447987003392, -1.169422842065824, -0.6657737930180916, 
0.6657736844488363,1.1694228046197082,1.3214479485830937}) 


The Positive Infinities of the main function will be find using the following method (tangent of the results 
): 


Tan[{-1.321447987003392, -1.169422842065824, -0.6657737930180916, 
0.6657736844488363,1.1694228046197082,1.3214479485830937}] 


Result : 


{-3.926991132598397, -2.3561946032087975, -0.7853982329054199, 
0.7853980573651929,2.356194357874862,3.9269905016888096} 



—oo (Negative Infinities) of the function can be find with the Initial guesses of 
{—1.32, —1.165, —.6,0.7,1.17,1.322} through the transformed function . 

The Mathematica code to find the Negative Infinities of the transformed function is as follows : 

FindRoot[7r/2 + ArcTan[Tan[2 * Tan[x]] 

+ Tan[x] A 3], {x, (-1.32, -1.165, -.6, 0.7, 1.17, 1.322}}, Maxlterations -» 300] 


Result : 


(x -» {-1.3214479677754711, -1.1694228247122458, -0.6657737474122314, 
0.6657737517571516,1.169422824845148,1.3214479678320301}} 


The Negative Infinities of the main function will be find using the following method (tangent of the 
results ): 


Tan[{-1.3214479677754711, -1.1694228247122458, -0.6657737474122314, 
0.6657737517571516,1.169422824845148,1.3214479678320301}] 


Result : 


{-3.9269908168517462, -2.3561944895141798, -0.7853981591675699, 
0.7853981661926553,2.356194490384908,3.9269908177805157} 


Now we want to go further and find the infinities that are not visible on the distorted plot with the 
domain of [— ~] • To find the invisible asymptotes we should use shorter domains for the infinity plot . 
We just show one example of invisible asymptote. So we plot the transformed function in a shorter 
domain as follow : 


Plot[ArcTan[Tan[2 * Tan[x]] + Tan[x] A 3], {x, —1.391, —1.3908}, PlotRange -> All] 


Result : 



+oo (Positive Infinities) of the function can be find with the Initial guesses of {—1.3909} through the 
transformed function . 

The Mathematica code to find the Positive Infinitiy of the transformed function is as follows : 

FindRoot[— 7t/2 + ArcTan[Tan[2 *Tan[x]] + Tan[x] A 3],{x, — 1.3909), Maxlterations 
— > 1000, AccuracyGoal -> 6] 

The Result is : 

(x -» -1.3908719929660696} 


The Positive Infinitiy of the main function will be find using the following method (tangent of the results 

): 

Tan[— 1.3908719929660696] 


The Result : 


-5.49778729840062 


—oo (Negative Infinities) of the function can be find with the Initial guesses of {—1.3909} through the 
transformed function . 

The Mathematica code to find the Negative Infinitiy of the transformed function is as follows : 


FindRoot[+7r/2 + ArcTan[Tan[2 *Tan[x]] + Tan[x] A 3],{x, — 1.3985}, Maxlterations 
— > 1000, AccuracyGoal -> 6] 

The Result is : 


x -» -1.390871972974362} 



The Positive Infinitiy of the main function will be find using the following method (tangent of the results 

): 


Tan[— 1.390871972974362] 


The Result : 


-5.49778667414632 





Developing A Method To Find The Absolute Maximum And 
Minimum Of Functions In The Domain (— 00 , + 00 ) . 


T o find the Absolute maximum and minimum of functions in the domain (—00, +00) we 
can use the ArcTan transformation and find the Absolute maximum and minimum of the 

transformed function in the domain + ~] • 


We introduce the method by some examples : 


Example 1: 


Find the Absolute Maximum and Minimum of the following function : 
f(x) = x 2 

The Maple Code to define the function and its transformed mode is as follows : 


fl •= X > x l 


f2 ■= x— >arctan(/7 (tan(x) ) ) 


x— 


x —> arctan (fl ( tan(x ) ) ) 

The Maple code to calculate the Maximum of the transformed function is as follows 


maximize 


f2(x),x=--^- location 


1 r 

rr 11 

1 



r 1 1 

1 

— n. 

S' 

11 

1 

■i 





, — 71 

2 ’1 

[[ 2 y 

2 J 



{ 2 \ 

1 2 


The Mathematica code to calculate the Absolute Maximum of the main function is as follows 


Tan § 


The result is : 


Complexlnfinity 



92 


The Maple code to check the result is as follows : 


maximize(fl (x), location ’ 


{ [ {x = co }, oo ], [ {x = - oo }, oo ] } 


As you Can see the result is true . 

Now we Calculate the absolute minimum : 

The Maple code to calculate the Minimum of the transformed function is as follows : 


minimize 


f2(x),x =--^- location J 


0,{[{x = 0},0]} 

The Maple code to calculate the Absolute minimum of the main function is as follows 


tan(O) 


0 


The Maple code to check the result is as follows : 

minimize(fl (x), location ) 

0, {[ {x = 0}, 0]} 

As you can see the result is true . 


Example 2 : 

Find the Absolute minimum and maximum of the following function : 

2 1 

v “ 


fix) = 


e — - 


The Maple Code to define the function and its transformed mode is as follows . 


fl ■= x- 


-x 2 1 

C 2 


e - * 2 L 

6 2 


f2 ■= x— >arctan(//(tan(x))) 


x— >arctan(/7 (tan(x) ) ) 


The Maple code to calculate the Maximum of the transformed function is as follows : 


maximize 


f2{x),x - location 


arctan I — |, 


{x = 0}, arctan^ yj J 


The Maple code to calculate the Absolute Maximum of the main function is as follows 

■(*)) 


tan arctan 


The Maple code to check the result is as follows : 

maximize(fl (x), location ) 


2 ’ 


{x = 0}, 


1 


As you Can see the result is true . 

Now we Calculate the absolute minimum : 

The Maple code to calculate the Minimum of the transformed function is as follows : 


minimize 


f2 (x),x =-~zr location 
J v ' 2 2 


0, { [ {x = - arctan (V ln(2) ) }, o], [ {x = arctan (V ln(2) ) }, o] } 

The Maple code to calculate the Absolute minimum of the main function is as follows : 

tan(0) 

0 

The Maple code to check the result is as follows : 


minimize{fl (x), location \ 


0,{[{* = Vta(2 n, o], [ {x = -JH2T], o]} 


As you can see the result is true . 
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Analitycal method 


N 


ow after seeing the applications of Infinity Plots , we will develop an Analytical method 
to find the unsolved limits of mentioned functions at the introduction of the article . 


The method is as follows : 


We will use the transformed function using the following transformation : 


ArcTan (/(rcm(x))) 


Df = 




n n 
2 ’ 2 . 


And we can use x — * (q — $MachineEpsilon) instead of the x — * +oo in the transformed 
function . 


So the limit will be changed to 


n lim [ ArcTan(f(Tan(x)yj] 

$MachineEpsilon) ^ ' 


Where/is the defined function : 


f=Fu action [{x},lntg[x]/x] 


So the Mathematica Code to Solve the Limit will be as follows : 


lntg=Function[{x}, lf[lntegerPart[x]>0,lntegerPart[x],lntegerPart[x]-l]]; 
F2= Function[{x},ArcTan[Fl[Tan[x]]]]; 

Y=F2[ j — $MachineEpsilon]; 

Tan[Y] 


Output : 


1 


So you can see that the Analytical method and the graphical method yield the same result . 


3d infinity plots 


/ 


t is possible to develop 3D Infinity Plots . The technique is the same as the 2D Infinity Plots . 
Here is some examples : 


F= Function[{x,y},{x A 2+y A 2,-x A 2-y A 2}]; 

Plot3D[ArcTan[F [Tan[x],Tan[y]]],{x,-7r /2, n /2},{y,- n /2, n /2},ColorFunction-»"RustTones"] 



F= Function[{x,y},{x A 2-y A 2}] 

Plot3D[ArcTan[F [Tan[x],Tari[y]]],{x,- n /2, n /2} ; {y,- n /2, n 
/2},RegionFunction->Function[{x,y,z},2<x A 2+y A 2<9],BoxRatios->Automatic] 


1 



F= Function[{x,y},x A 2-y A 2] 

F5Region= Function[{x,y,z},2<x A 2+y A 2<9] 

ContourPlot[ArcTan[F5[Tan[x],Tan[y]]],{x,- n /2, n /2},{y,- n /2, n 
/2},RegionFunction->Function[{x,y,z},2<x A 2+y A 2<9]] 





System of nonlinear equations 



/ n this part of the article we are going to develop some code to solve the system of 
nonlinear equations through extended Broyden method . The Broyden LUA codes are 
copyrighted by PAP and the extended Broyden method codes are copyrighted by the author 
of this article . 

One of my friends with the nickname of Mohammad suggested me to extend the broyden 
method and use the Infinity applications that we discussed earlier in the method to have an 
improved method to solve the system of non-linear equations . After thinking about the 
method I could find an algorithm to find the solutions of the system of non-linear equations. 
First we will describe the Broyden method then we will extend it . 


• Broyden/^ 

The function Broyden solves a system of non-linear equations by implementing Broyden's 
method . This method is iterative, and tries to find a solution, starting by a user-supplied initial 
guess . It is often considered as the method of choice for solving systems of non-linear 
equations . 

Syntax/ 2 ^ 


Root=Broyden(xguess, F,J, eps,maxit,show) 


returns a vector root,which is the solution of a system of non-linear equations . The general 
form of a system of N non-linear equations with N unknowns x lt x 2 , ... , x N is F(x)=0 , which is 
expressed analytically as 

Fi(x 1 ,x 2 , ■■■ , %n) = 0 
F 2 Oi,* 2 , ...,x n ) = 0 


Fn(. x i> x 2 , ...,Xn) — 0 


Xguess is initial guess of the solution , and should be a vector of N elements . F is the name of a 
user-defined function F(x) , which returns a vector defining the left-hand-side of the equations 
to be solved , given a vector x defining the independent variables . The arguments J,eps,maxit , 
and show are optional . 

The argument J must be the name of a function J(x) defining the Jacobian matrix J(x) . If omitted 
, this argument takes the default value J="auto" , which means that Jacobian will be computed 
numerically . 



Eps defines the desired accuracy ( default : Epsilon , i.e , 1. 12 x 10 -16 ); maxit defines the 
maximum number of iterations (default : 20) ; show is a Boolean argument that controls 
whether progress of calculation will be displayed or not ( default : false ) . 

Example^ 

Consider the system of non-linear equations 

X 1 ~ x 2 + X 3 ~ 5 = 0 

x\ + x 2 - 3 = 0 

x 1 + x 2 x 3 = 0 


In this case , the multivariate function F(x) is defined as 


Hx) = Hxi,x 2 ,x 3 ) = 


*1 -x 2 + x 3 - 
x \ + x 2 - 3 
x ± + x 2 x 3 



Its Jacobian is 


1 


J(x) = i(Xi, x 2 , x 3 ) = 


2x ± 


1 


-4x1 1 

3x\ 0 

x 3 2x ± x 2 . 


The example program XBroyden uses the function Broyden to find the solution of this system 
of non-linear equations with initial guess x t = 2 , x 2 = 0 , x 3 = 3 . In order to check the result 
, the program computes the function values at the solution obtained . 


require("LNA/Broyden ") 
local function F(x) 

return {x[l ]-x[2] A 4+x[3]-5,x[l ] A 2+x[2] A 3-3,x[l ]+x[2]*x[3] A 2} 
end 

local function J(x) 

return {{1, -4 *x[2] A 3, 1 }, {2 *x[l ],3*x[2] A 2,0}, {l,x[3 ] A 2,2*x[2 ] *x[3]}} 
end 

local root 

root=Broyden({2, 0, 3 }, F,J) 
print("Solution:") 

for i,v in ipairs(root) do print(i,v) end 

print("\nFunction value:") 

for i,v in ipairs(F(root)) do print(i,v) end 


Example Program : XBroyden 


The Figure bellow shows the results obtained by running this program . Note that , in this 
example , the Jacobian is defined by the user . Alternativley , one can call the function Broyden 
as root=Broyden({2,0,3},F) . In this case , the Jacobian is omitted , and it is computed 
numerically. 

The reader should modify the initial guess in the above example , to see how the solution 
obtained is affected by a "bad" initial guess . For example , try following cases : 

1. xguess={3,l,6} is not a good initial guess , but the solution obtained is correct . 

2. xguess={l,0,5} is not a good initial guess , but the solution obtained is accurate , 
although the maximum number of iteration is reached , and the accuracy criterian is not 
satisfied . 

3. xguess={0,l,0} is a very bad initial guess ; the maximum number of iteration is reached , 
the accuracy criterion is not satisfied ,and the "solution" obtained is wrong . 

Remarks/ 2 ^ 

A good choice of initial guess xguess is crucial for Broyden's algorithm . Although you may get 
correct results using an initial guess far from the correct solution , it is always better to supply a 
good initial guess . Failure to do so may cause Broyden to return a totally wrong "solution" ; in 
this case , a nil value will be returned , and you should alter the initial guess. 

Usually , there is no need to change the default value for "accuracy" eps . The algorithm should 
converge after a few iterations ; if i t does not , it is almost sur e that it will never converge , 

V Console Break 

Solution: 

1 1 . 7332857423775 

2 -0. 16235405884328 

3 3. 267403046091 

Function value: 

1 0 

2 4. 440S920935006e-016 

3 0 
Done. 

Press [<-] or [EXE]... 


Done. ijm] 

Results obtained by the example program XBroyden . 

Due to wrong initial guess , or wrong input data . In other words , you will probably never need 
to use the optional argument maxit . The most common error when using Broyden is to define 
F incorrectly ; is most cases , this will cause wrong results , obtained after the maximum 
number of iterations is reached . 

If the Jacobian is omitted , it is computed numerically by using the function Jacobian . Although 
numerical computation of the Jacobian is usually accurate enough , it is better to supply it by 
using the argument . In most cases , computing the Jacobian analytically is not difficult at all . 



Extended Broyden Method 


Now we are going to develop the extended Broyden method to solve the system of non-linear 
equations to find the solutions of the system as the variables are varied in the Range (—00, +co). 

The Broyden method uses initial guesses . we can improve the method by increasing and decreasing 
the initial guesses in a domain of a transformed function. We will introduce the method through an 
example . 

Assume that we have the following function ( it is the example that we used for introducing the 
Broyden method in the previous section ): 

local function F(x) 

return {x[l ]-x[2] A 4+x[3]-5,x[l ] A 2+x[2] A 3-3,x[l ]+x[2]*x[3] A 2} 

end 


First we should transform the main function using the following Formula : 


ArcTan(F(Tan(x)) 


Then we should increase and decrease the variables that are included in the initial guesses just in the 
domain of ■ 

The transformed function in lua code will be as follows : 


local function FTrans(x) 

return math.atan(F(math.tan(x))) 

end 


Instead of using ArcTan in the above code , we will use a power series representing the ArcTan function 
.The Mathematica code to represent ArcTan as a power series with good accuracy is as follows : 

Normal[Series[ArcTan[x], {x, 0, 50}]] 


The result will be: 


r 3 r 5 v 7 v 9 r ll r 13 r 15 r 17 v 19 r 21 r 23 r 25 

A A A A A A A A A A A A 

F = Fun«io„[{*},* — T + T - y + — + — - — + — 

y.21 y.29 V 31 v 33 v 35 ^37 -y-39 V 41 <y43 v 45 y 47 v 49 

A A A A A A A A A A A A 

— "27" "29 3T "33 35" "37 39" "41" — "43" "45" — "47" 49"^ 


So the terms of Ftrans function will be as follows : 



Term 1 : 

The mathematica Code to produce Ftrans term 1 is : 
F[Tan[x[l]] - (Tan[x[2]]) A 4 + Tan[x[3]] - 5] 


the result will be : 


-5 + Tan[x[l]] - Tan[*[2]] 4 + Tan[x[3]] - - (-5 + Tan[x[l]] - Tan[x[2]] 4 + Tan[x[3]]) 3 

1 1 
+ - (-5 + Tan[jc[l]] - Tan[x[2]] 4 + Tan[x[3]]) 5 - - (-5 + Tan[x[l]] - Tan[x[2]] 4 + Tan[x[3]]) 7 

1 1 
+ - (-5 + Tan[x[l]] - Tan[x[2]] 4 + Tan[x[3]]) 9 - — (-5 + Tan[x[l]] - Tan[x[2]] 4 + Tan^]]) 1 


+ (— 5 + Tan[x[l]] 

+ W (_5 + Tan[x[1]] 

+ -(-5 + Tan[x[l]] 
+ ^(-5 + Tan[*[l]] 
+ -(-5 + Tan[x[l]] 

+ ^ (_5 + Tan[x[1]] 
+ ^(-5 + Tan[x[l]] 

+ ? 1 (_5 + Tan[x[1]] 
+ ^ ( -5 + Tan[ * [1]] 


11 

- Tan[x[2]] 4 + Tan[x[3]]) 13 - (-5 + Tan[x[l]] - Tan[x[2]] 4 + Tan[x[3]]) 15 

- Tan[x[2]] 4 + Tan[x[3]]) 17 --^(-5 + Tan[x[l]] - Tan[x[2]] 4 + Tan[x[3]]) 19 

- Tan[x[2]] 4 + Tan[x[3]]) 21 - -^(-5 + Tan[x[l]] - Tan[*[2]] 4 + Tan[*[3]]) 23 

- Tan[x[2]] 4 + Tan[x[3]]) 25 - -^(-5 + Tan[x[l]] - Tan[*[2]] 4 + Tan[*[3]]) 27 

- Tan[x[2]] 4 + Tan[x[3]]) 29 - ^-(-5 + Tan[x[l]] - Tan[x[2]] 4 + Tan[x[3]]) 31 

- Tan[x[2]] 4 + Tan[x[3]]) 33 - -^(-5 + Tan[x[l]] - Tan[x[2]] 4 + Tan[x[3]]) 35 

- Tan[x[2]] 4 + Tan[x[3]]) 37 - -^(-5 + Tan[x[l]] - Tan[x[2]] 4 + Tan[x[3]]) 39 

- Tan[x[2]] 4 + Tan[x[3]]) 41 - i(-5 + Tan[x[l]] - Tan[x[2]] 4 + Tan[x[3]]) 43 

- Tan[x[2]] 4 + Tan[x[3]]) 45 -. — (-5 + Tan[x[l]] - Tan[x[2]] 4 + Tan[x[3]]) 47 


+ — (-5 + Tan[x[l]] - Tan[x[2]] 4 + Tan[x[3]]) 4 


Term 2 : 

The mathematica Code to produce Ftrans term 2 is : 
F[(Tan[x[l]]) A 2 - (Tan[x[2]]) A 3 - 3] 



the result will be : 


— 3 + Tan[*[l]] 2 — Tan[x[2]] 3 — — (— 3 + Tan[x[l]] 2 — Tan[x[2]] 3 ) 3 + — (— 3 + Tan[x[l]] 2 — Tan[x[2]] 3 ) 5 

- ^ (-3 + Tan[x[l]] 2 - Tan[x[2]] 3 ) 7 + ^ (-3 + Tan[x[l]] 2 - Tan[x[2]] 3 ) 9 

- jj(— 3 + Tan[x[l]] 2 - Tan[x[2]] 3 ) n + -^(-3 + Tan[x[l]] 2 - Tan[x[2]] 3 ) 13 

- (-3 + Tan[x[l]] 2 - Tan[x[2]] 3 ) 15 + (-3 + Tan[x[l]] 2 - Tan[x[2]] 3 ) 17 

- ^ ( -3 + Tan[x[l]] 2 - Tan[x[2]] 3 ) 19 + -^(-3 + Tan[*[l]] 2 - Tan[*[2]] 3 ) 21 

- ^ ( -3 + Tan[x[l]] 2 - Tan[x[2]] 3 ) 23 + ^ (-3 + Tan[x[l]] 2 - Tan[*[2]] 3 ) 25 

- ^ ( -3 + Tan[x[l]] 2 - Tan[x[2]] 3 ) 27 + ^ (-3 + Tan[x[l]] 2 - Tan[*[2]] 3 ) 29 

- (-3 + Tan[x[l]] 2 - Tan[x[2]] 3 ) 31 + (-3 + Tan[x[l]] 2 - Tan[x[2]] 3 ) 33 

- (-3 + Tan[x[l]] 2 - Tan[x[2]] 3 ) 35 + (-3 + Tan[x[l]] 2 - Tan[x[2]] 3 ) 37 

“ (-3 + Tan[x[l]] 2 - Tan[x[2]] 3 ) 39 + ^ (-3 + Tan[x[l]] 2 - Tan[x[2]] 3 ) 41 

- — (—3 + Tan[x[l]] 2 - Tan[x[2]] 3 ) 43 + — (—3 + Tan[x[l]] 2 - Tan[x[2]] 3 ) 45 

- i (-3 + Tan[x[l]] 2 - Tan[x[2]] 3 ) 47 + ^ (-3 + Tan[x[l]] 2 - Tan[*[2]] 3 ) 49 


Term 3 : 

The mathematics Code to produce Ftrans term 3 is : 
F[(Tan[x[l]]) A 2 - (Tan[x[2]]) A 3 - 3] 


the result will be : 

Tan[x[l]] + Tan[x[2]]Tan[x[3]] 2 — -(Tan[x[l]] + Tan[x[2]]Tan[;e[3]] 2 ) 3 + -(Tan[x[l]] + Tan[x[2]]Tan[;e[3]] 2 ) 5 

1 1 

— -(Tan[x[l]] + Tan[x[2]]Tan[x[3]] 2 ) 7 + -(Tan[x[l]] + Tan[x[2]]Tan[x[3]] 2 ) 9 

- (Tan[x[l]] + Tan^pJJTan^]] 2 ) 11 + i (Tan[x[l]] + Tan[x[2]]Tan[x[3]] 2 ) 13 

- (Tan[x[l]] + Tan[x[2]]Tan[*[3]] 2 ) 15 + (Tan[x[l]] + Tan[x[2]]Tan[*[3]] 2 ) 17 

— ^j(Tan[x[l]] + Tan[x[2]]Tan[*[3]] 2 ) 19 + ^-(Tan[x[l]] + Tan[*[2]]Tan[x[3]] 2 ) 21 
“^(Tan^I 1 ]] + Tan[x[2]]Tan[x[3]] 2 ) 23 + ^(Tan[*[l]] + Tan[x[2]]Tan[x[3]] 2 ) 25 
-^(Tanlxfi]] + T a n [x[2]] T an[x[3]] 2 ) 27 + ^('ran[x[ 1 ]] + T an[x[2]]Ta n [x[3]] 2 ) 29 
-^-(Tan[x[l]] + Tan[x[2]]Tan[x[3]] 2 ) 31 +^(Tan[x[l]] + Tan[x[2]]Tan[x[3]] 2 ) 33 

- (Tan[x[l]] + Tan[x[2]]Tan[x[3]] 2 ) 35 + ^(Tan[x[l]] + Tan[x[2]]Tan[x[3]] 2 ) 37 

— — (Tan[x[l]] + Tan[x[2]]Tan[x[3]] 2 ) 39 + — (Tan[x[l]] + Tan[x[2]]Tan[x[3]] 2 ) 41 
~^(Tan[x[l]] + Tan[x[2]]Tan[x[3]] 2 ) 43 + ^(Tan[x[l]] + Tan[x[2]]Tan[x[3]] 2 ) 45 
--^(Tan[*[l]] + Tan[x[2]]Tan[x[3]] 2 ) 47 + -^(Tan[*[l]] + Tan[x[2]]Tan[x[3]] 2 ) 49 


In this example we have 3 variables and each variable can increase or decrease in the domain of 



So we will have 3 loops that in the first loop X[l] variable will be changed in the domain with the 

user defined step (accuracy) and X[2] and X[3] will be — | and in the second loop X[2] will be changed in 
the domain [ — while X[3] is — ^and X[l] will be changed in the domain for each X[2] value 

and in the last loop X[3] will be changed in the domain while X[2] will be changed in the domain 

f ’f] ^ or eac ^ ^[3] value and X[l] will be changed in the domain [ — f° r eac h X[2] value . 

In each loop we will see if the try will converge to a solution with satisfied accuracy or not . If the 
solution is converged with satisfied accuracy we will use the result as a good initial guess for the 
Broyden method with the main function . when we reach to the end of first loop the next loop will be 
started and we will continue upto the end of the procedure . 

In this method the solutions for the Ftrans function will be the initial guesses for the F function . 

Using this method the real solutions of the system of equations with each variable varied in the domain 
(—oo, +oo) w j|| be optained . 

Syntax 

BroyINFfnvar, F, Ftrans, step, J,Jtrans,eps,maxit, show) 

Returns a vector root , which is the solution of a system of non-linear equations . The general 
form of a system of N non-linear equations with N unknowns x lt x 2 , ... , x N is F(x)=0 , which is 
expressed analytically as 


F iOi,*2, -,X N ) = 0 
F 2 (.Xi,X 2 , ■■■ , Xf^) = 0 


Fn(x\,x 2 , ... ,x N ) = 0 

nvar is the maximum number of variables is the equations . F is the name of a user-defined function F(x) 
, which returns a vector defining the left hand side of the equations to be solved , given a vector x 
defining the independent variables . Ftrans is the transformed function of F using the following 
transformation . 


ArcTan(F(Tan(x))) 


step is the step size in the solution , can be started from 0.01 to the Epsilon . smaller step size take 
more time to solve the system of equations and more results may be found . 

The arguments J , Jtrans , eps , maxit and show are optional . 

The argument J must be the name of a function J(x) defining the Jacobian matrix J(x) . If omitted 
, this argument takes the default value J="auto" , which means that Jacobian will be computed 
numerically . 

The argument Jtrans must be the name of a function Jtrans(x) defining the Jacobian matrix 
Jtrans(x) . If omitted , this argument takes the default value Jtrans="auto" , which means that 
Jacobian of the transformed function will be computed numerically . 

Eps defines the desired accuracy ( default : Epsilon , i.e , 1. 12 x 10 -16 ); maxit defines the 
maximum number of iterations (default : 20) ; show is a Boolean argument that controls 
whether progress of calculation will be displayed or not ( default : false ) . 

Example Driver program 

The Driver program to solve the example will be as follows : 


require("LNA/BroylNF") 


local function F(x) 

return {x[l ]-x[2] A 4+x[3]-5,x[l ] A 2+x[2] A 3-3,x[l ]+x[2]*x[3] A 2} 
end 


local function Ftransformer(x) 
return 

x-x A 3/3+x A 5/5-x A 7/7+x A 9/9-x A ll/ll+x A 13/13-x A 15/15+x A 17/17-x A 19/19+x A 21/21- 
x A 23/23+x A 25/25-x A 2 7/27+x A 29/29-x A 31/31+x A 33/33-x A 35/35+x A 3 7/3 7-x A 39/39+x A 41/41- 
x A 43/43+x A 45/45-x A 4 7/4 7+x A 49/49 
end 

local function Ftrans(x) 
return 

{Ftransformer(F({math. tan(x[l] ), math. tan(x[2]), math. tan(x[3] )})[!]), Ftransformer(F({math.tan(x[ 
1]), math. tan(x[2]), math. tan(x[3])}) [2] ),Ftransformer(F({math.tan(x[l]), math. tan(x[2]), math. tan( 
x[3])})[3])} 
end 


BroyINFf 3, F, Ftrans, 0.3) 


The results will be as follows : 


Console Break 


Console Break 


# Initial Guesses = 1331 £ 

k 

Solution 2 : 

1 -1.7320346936097 

2 0.033217613493455 

3 6.7320363269222 


2 Result(s> Found With 

The Defined Accuracy . 

They Are As Follows : 


I 

Solution 1 : 

1 1 . 7332357423775 

2 -0. 16235405334323 

3 3. 267409046091 

Function value: 

1 0 

2 4.4403320385006e- 
016 

3 0 

g 


Function value: 

1 0 

2 0 

3 0 

Done. 



r 

Press [<-] or [EXE]... 

T 

Done. ijm] 


Done. ijm] 



The Extended Broyden method Code and its requirements : 


BroyINF 


requiref "table ", "LNAutils/Epsilon ", "LNAutils/Pi", ", LNAutils/MatMul ", "LNAutils/MatTrans", "LNA/J 
acobia ", "LNA/LUdecom ", "LNA/Broyden ") 

function BroylNF(nvar, F, Ftrans,step, ...) 
local J=" auto" 
local Jtrans-'auto" 
local eps=Epsilon 
local maxit=20 
local show=false 
if arg["n"]>=l then 
J=arg[l] 

if arn["n "]>=2 then 
Jtrans=arg[2] 
if arg["n"]>=3 then 
eps=arg[3] 
if arg["n"]>=4 then 
maxit=arg[4] 

ifarg["n"]>=5 then 



show=arg[5] 

end 

end 

end 

end 

end 

local n=nvar 
local stepi=l 

local xguessvalue = {} 
if (n == 1) then 
for all=-Pi/2, Pi/2, step do 
xguessvalue[stepi] = {all} 
stepi=stepi+l 
end 
end 


if (n == 2) then 
for al2=-Pi/2, Pi/2, step do 
for all=-Pi/2, Pi/2, step do 

xguessvalue[stepi] = {all, al2} 
stepi=stepi+l 
end 
end 
end 


if (n == 3) then 
for al3=-Pi/2, Pi/2, step do 
for al2=-Pi/2, Pi/2, step do 
for all=-Pi/2, Pi/2, step do 
xguessvalue[stepi] = {all, al2, al3 } 
stepi=stepi+l 
end 
end 
end 
end 


if (n == 4) then 
for al4=-Pi/2, Pi/2, step do 
for al3=-Pi/2, Pi/2, step do 
for al2=-Pi/2, Pi/2, step do 
for all=-Pi/2, Pi/2, step do 


xguessvalue[stepi] = {all, al2, al3, al4} 
stepi=stepi+l 

end 

end 

end 

end 

end 

if (n==5) then 
for al5=-Pi/2, Pi/2, step do 
for al4=-Pi/2, Pi/2, step do 
for al3=-Pi/2, Pi/2, step do 
for al2=-Pi/2, Pi/2, step do 
for all=-Pi/2, Pi/2, step do 

xguessvalue[stepi] = {all, al2, al3, al4} 

stepi=stepi+l 

end 

end 

end 

end 

end 

end 


lengthguess=stepi-l 

print/'# Initial Guesses = "..lengthguess) 

local root={} 

local solution={} 

local solutionc={} 

for ii = 1, lengthguess do 

solution [ii]=Broyden(xguessvalue[ii],Ftrans,Jtrans) 
if (solution[ii] ~=nil) then 
solutionc = {}; 
for i=l,n do 

solutionc[i]=math.tan(solution[ii][i]); 

end 

root[ii]=Broyden(solutionc,F,J,eps,maxit,show) 

else 

root[ii]=nil 

end 

end 



result = {}; 
result i = 0; 
mresult = {}; 
stepi2last = stepi-1; 

for steppass2i=l,stepi2last do 
for m=l,(2 A n) do 

if root[steppass2i]~=nil then 
pvalue = 0; 

for steppass3i=l,steppass2i do 
if (steppass3i ~= steppass2i ) then 
mprocess = (2 A n) 

else 

mprocess = m; 

end 

for m2=l, mprocess do 
pkey = 0; 

if root[steppass3i]~=nil then 
for i=l,n do 

if(math.abs(root[steppass2i][i] - root[steppass3i][i]) < 0.0000000001) then 
pkey = pkey + 1 

end 

end 

end 

if (pkey == n) then 
pvalue = pvalue + 1 

end 

end 

end 

if (pvalue == 1) then 
result i = result i + 1 
result [ resulti ]=steppass2i 
mresult[resulti]=m 

end 

end 

end 

end 


p r/n tr- ===================== 'V 

if resulti==0 then 

printf'No Result Found With The Defined Accuracy") 
else 



print(resulti.. " Result(s) Found With The Defined Accuracy . They Are As Follows :") 

print(" ") 

for istep=l, resulti do 
print("Solution "..istep.." :") 
for i,v in ipairs(root[result[istep]]) do print(i,v) end 
print("\nFunction value:") 

for i,v in ipairs(F(root[result[istep]])) do print(i,v) end 

print!"- ") 

end 

end 


end 

export{BroylNF=BroylNF} 


Broyden 


require("table", "LNAutils/Epsilon", "LNAutils/MatMul", "LNAutils/MatTrans", "LNA/Jacobia", "LNA 
/LUdecom") 

function Broyden(xguess,F,...) 
local J=" auto" 
local eps=Epsilon 
local maxit=20 
local show= false 
if arg["n"]>=l then 
J=arg[l] 

if arg["n"]>=2 then 
eps=arg[2] 
if arg["n"]>=3 then 
maxit=arg[3] 
if arg["n"]>=4 then 
show=arg[4] 
end 
end 
end 
end 

-Numerical Jacobian: 
if J==" auto" then 
local function Jnum(x) 
return Jacobian!F,x ) 



end 

J=Jnum 

end 

-Local variables: 
local n=#xguess 

local LU,indx,v,vmax,Fcur,Fpre,Ainv,Ainvy,DTAinv,DTAinvy,DmAinvy,aux 

local y={} 

local DmAinvy={} 

local x=table.copy(xguess) 

-First root estimation: 

Fpre=F(x) 

LU,indx=LUdecompose(J(x)) 

Ain v=L Uin verse (L U, indx) 
v=MatMul(Ainv,Fpre) 
for i=l,n do 
x[i]=x[i]-v[i] 
end 

if show then 

print('' i max discrepancy") 
vmax=math. abs( v[l ]) 
for i=2,n do 
aux=math. abs( v[i]) 
if aux>vmax then 
vmax=aux end 
end 

printf("%2i %e\n",l,vmax ) 
end 

for iter=2,maxit do 
-Jacobian inverse: 

Fcur=F(x) 
for i=l,n do 
y[i]=Fcur[i]-Fpre[i ] 
end 

Ainvy=MatMul(Ainv,y) 

DTAinv=MatMul(v / Ainv) 
for i=l,n do 

Dm A in vy[i]=-v[i]-A in vy[ i ] 
end 

aux=MatMul(MatTrans(DmAinvy),{DTAinv}) 

DTAinvy=MatMul(DTAinv / y) 
for i=l,n do 
for j=l,n do 

Ainv[ i ][j]=Ainv[ i ][j ]+aux[ i ][j ]/DTAinvy 
end 
end 

-Next root estimation: 




v=MatMul(Ainv,Fcur) 
for i=l,n do 
x[i]=x[i]-v[i] 
end 

-Loop or return: 
vmax=math. abs( v[l ]) 
for i=2,n do 
aux=math. abs( v[i]) 
if aux>vmax then 
vmax=aux end 
end 

if show then 

printf("%2i %e \n ",iter, \ /max) 
end 

if vmax<eps then 
return x 
else 

Fpre=Fcur 

end 

end 

-Return with warning: 

return nil 
end 

export{Broyden=Broyden } 


Epsilon 


Epsilon=1.12E-16 

export{Epsilon=Epsilon} 


LNAutils/MatMul 


require(''LNAutils/MatCol", "LNAutils/MatTrans ") 
function MatMul(A,B) 
local C={} 

local rowA,colA,rowB,colB 
if type(A[l])=="number" then 
A={A} 




end 

rowB=#B 

if type(B[l])=="number" then 
B=MatTrans(B);colB=l 
else 

colB=#B[l] 

end 

rowA=#A;colA=#A[l] 
if colA~=rowB then 

print("Error: MatMul: Undefined matrix multiplication. ") 
return nil 
end 

for i=l,rowA do 

CUM 

for j=l,colB do 

cim=o 

for k=l,colA do 
C[i][j]=C[i][j]+A[i][k]*B[l<][j] 
end 
end 
end 

if colB==l then 
C=MatCol(C,l) 
end 

if row A==l then 
C=C[1] 
end 

return C 
end 

export{MatMul=MatMu 1} 


LNAutils/MatCol 


function MatCol(A,col) 
local C={} 

iftype(A[l])=="table” then 
for row=l,#A do 
C[row]=A[row][col] 
end 
else 

C[l]=A[col] 



end 

return C 
end 

export{MatCol=MatCol} 


LNAutils/MatTrans 


require(”LNAutils/MatCol") 
function MatTrans(M) 
local Mt={} 
local cols 

if type(M[l])==''table” then cols=#M[l] 

else cols=#M end 

for i=l,cols do 

Mt[i]=MatCol(M,i) 

end 

return Mt 
end 

export{MatTrans=MatTrans} 


LNA/Jacobia 


require("table") 
function Jacobian(F,x) 
local n=#x 
local Fx=F(x) 
local dx={} 
local Fxk={} 
local xaux 
local J={} 
for k=l,n do 
ifx[k]~=0 then 
dx[k]=lE-8*x[k] 
else 

dx[k]=lE-8 

end 

xaux=table. copy(x) 
xaux[k]=x[k]+dx[k] 


Fxk[k]=F(xaux) 

end 

for i=l,n do 

m={. i 

for k=l,n do 

J[i][k]=(Fxk[k][i]-Fx[i])/dx[k] 

end 

end 

return J 
end 

export{Jacobian=Jacobian } 


LNA/ LUdecom 


requiref "table ", "LNAutils/MatCol", "LNAutils/Matldent", "LNAutils/MatTrans ") 
function LUdecompose(A) 
local LU={} 
local indx={} 
local parity=l 
local n=table.getn(A) 
local p,Dp,Dc,aux 
for i=l,n do 
L U[i]=table. copy(A[i ] ) 
indx[i]=i 
end 

for j=l,n-l do 
-Upper triangular: 
for i=l,j do 
for k=l,i-l do 

LU[i][j]=LU[i][j]-LU[i][k]*LU[k][j] 

end 

end 

-Lower triangular: 
for i=j+l,n do 
for k=l,j-l do 

LU[i][j]=LU[i][j]-LU[i][k]*LU[k][j] 

end 

end 

-Pivoting: 

P=j 

Dp=math.abs(LU[p][j]) 
for i=j+l,n do 


Dc=math.abs(LU[i][j]) 
if Dc> Dp then 
p=i;Dp=Dc 
end 
end 

if Dp==0 then 

print(''Warning: LUdecompose: Singular matrix. ") 
return nil, nil, nil 
end 

-Permutation: 
if P~=j then 
for k=l,n do 

LU[j][k],LU[p][k]=LU[p][k],LU[j][k] 

end 

indx[j],indx[p]=indx[p],indx[j] 

parity=-parity 

end 

-Sub-diagonal division: 
uux=LU[j][j] 
for i=j+l,n do 
LU[i][j]=LU[i][j]/aux 
end 
end 

-Last column: 
for i=l,n do 
for k=l,i-l do 

L U[i][n ]=L U[i][n i ]-L U[i][k ] *L U[k][n ] 
end 
end 

return LU,indx, parity 
end 

function LUsubstitute(LU,indx,B) 
local X={} 
local Z={} 
local n,s 
if LU~=nil then 
n=table.getn(indx) 

-Lower system solution: 

Z[l]=B[indx[l]] 
for i=2,n do 
s=0 

for k=l,i-l do 
s=s+LU[i][k]*Z[k] 
end 

Z[i]=B[indx[i]]-s 


end 

-Upper system solution: 
X[n]=Z[n]/LU[n][n] 
for i=n-l, 1,-1 do 
s=0 

for k=i+l,n do 
s=s+LU[i][k]*X[k] 
end 

X[i]=(Z[i]-s)/LU[i][i] 

end 

return X 
else 

return nil 
end 
end 

function LUdeterminantfLU, parity) 
local s=l 
if LU~=nil then 
for i=l, table. getn(LU) do 
s=s*LU[i][i] 
end 

return parity*s 
else 

return 0 
end 
end 

function LUinverse(LUJndx) 
local Ainv={ } 
local n,!d 
if LU~=ni I then 
n=table.getn(LU) 
ld=Matldent(n) 
for i=l,n do 

Ainv[i]=LUsubstitute(LU,indx,MatCol(ld,i)) 

end 

return MatTrans(Ainv) 
else 

return nil 
end 
end 

function LUsolve(A,B) 
local X,LU,indx, parity 
L U, indx / parity=LUdecompose(A) 

X=L Usubstitutef L U, indx, B ) 


return X 
end 

export{LUsolve=LUsolve,LUdecompose=LUdecompose / LUsubstitute=LUsubstitute,LUdeterminant 
=L Udeterminant, L Uin verse=L Uin verse } 


LNAutils/Matldent 


function Matldent(n) 
local l={} 
for i=l,n do 

l[i]=0 

for j=l,n do 

nm=o 

end 

l[i][i]=l 

end 
return I 
end 

export{Matldent=Matldent} 


Another Example 

In this example we will have a function with 2 variable . The function is : 


local function F(x) 

return {x[l ]-x[2] A 2+5,x[l ]-math.sin(x[2])} 
end 


The Driver program 


require(”LNA/BroylNF") 

local function F(x) 

return {x[l ]-x[2] A 2+5,x[l ]-math.sin(x[2])} 
end 

local function Ftransformer(x) 
return 



x-x A 3/3+x A 5/5-x A 7/7+x A 9/9-x A ll/ll+x A 13/13-x A 15/15+x A 17/17-x A 19/19+x A 21/21- 
x A 23/23+x A 25/25-x A 2 7/27+x A 29/29-x A 31/31+x A 33/33-x A 35/35+x A 3 7/3 7-x A 39/39+x A 41/41- 
x A 43/43+x A 45/45-x A 47/47+x A 49/49 
end 


local function Ftrans(x) 
return 

{Ftransformerf F({math.tan(x[l]), math. tan(x[2] )})[l]),Ftransformer(F({math.tan(x[l]), math. tan(x[ 
2})})[2])} 
end 


BroyINFf 2, F, Ftrans, 0.1) 


Results 


V Console Break 


# Initial Guesses = 1024 


2 Result(s> Found With 
The Defined Accuracy . 
They Are As Follows : 


Solution 1 : 

1 -0.89851782355808 

2 -2.0252116374448 

Function value: 

1 0 
2 0 


Solution 2 : 

1 0.6866827734478 

2 2.3846766601466 

Function value: 

1 -8. 3317341970013e- 
016 

2 0 


Done. Gy] 


Checking the result through Mathematica : 


FindRoot[{y == (x) A 2 -5,y== Sin[x]}, {{x, 2}, {y, 1}}] 
FindRoot[{y == (x) A 2 - 5, y == Sin[x]}, {{x, -2}, {y, -1}}] 



Result : 


{x -> 2.38468, y -> 0.686683} 
{x -> -2.02521, y -> -0.898518} 


Plotting the functions using Mathematica: 


Plot[{ArcTan[(Tan[x]) A 2 - 5],ArcTan[Sin[Tan[x]]]}, {x, -([Pi]/2), [Pi]/2},PlotRange -> [Pi]/2] 


Result : 
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